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gaseous  absorption  within  the  LOWTRAN/ FAS CODE  formalism. 
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mation  to  numerical  multiple  scattering  treatments  for  solar,  thermal,  and 
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the  path  dependent  radiance  than  do  their  nonscattering  counterparts.  To 
Insure  consistency,  multiple  scattering  versions  of  LOWTRAN  and  FASCODE  are 
compared  to  one  another. 

Finally,  recommendations  are  presented  to  further  exploit  the  capabili¬ 
ties  introduced  here. 
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1 .  INTRODU CT ION 


I 


P 


1 . I  Background 

Optimal  design  and  deployment  of  electro-optical  (EO)  remote  sensing  and 
communication  systems  requires  accurate  modeling  and  prediction  of  the  effects 
of  the  ambient  environment  on  atmospheric  transmission.  Atmospheric  transmit¬ 
tance/radiance  models  such  as  AFGL’s  LOWTRAN  (Kneizys  et  al.,  1983)  and 
FASCODE  (Clough  et  al.,  1986)  have  been  developed  within  this  context  to  pro¬ 
vide  the  capability  to  assess  potential  adverse  environmental  impacts  on  EO 
system  performance. 

In  order  to  accurately  predict  atmospheric  effects  on  the  propagation  of 
visible,  infrared  and  microwave  radiation,  it  is  necessary  to  treat  the  ex¬ 
tinction  mechanisms  including  molecular  scattering  and  absorption,  and  particle 
(aerosols,  clouds  and  precipitation)  scattering  and  absorption,  characterizing 
the  ambient  atmosphere.  In  the  present  AFGL  transmittance/radiance,  models, 
these  processes  are  adequately  included  in  the  treatment  of  path  transmission. 
However,  simplified  treatments  are  employed  to  simulate  the  effects  of  scat¬ 
tering  on  the  calculation  of  radiance.  For  thermal  infrared  and  microwave 
radiation,  for  example,  particle  scattering  in  LOWTRAN  has  been  treated  as  an 
enhancement  to  extinction  but  not  as  a  source  term.  This  approach  leads  to  an 
underestimate  of  radiance  for  paths  where  multiple  scattering  is  important 
(Ben-Shalom  et  al.,  1980).  LOWTRAN  uses  the  single  scattering  approximation 
for  evaluating  solar  radiances  (Ridgway  et  al.,  1982),  While  the  single  scat¬ 
tering  implementation  is  straightforward,  its  application  introduces  errors 
which  are  functions  of  wavelength,  sun/sensor  geometry,  and  surface  optical 
properties  (cf,  Isaacs  and  Ozkaynak,  1980;  Dave,  1981).  These  errors  are  pri¬ 
marily  due  to  neglecting  higher  order  scattering  and  surface  reflection.  In 
FASCODE,  particle  scattering  has  been  treated  as  equivalent  to  absorption. 

All  scattered  radiation  is  thus  re-emitted  as  if  it  were  absorbed,  i.e.,  the 
scattered  radiation  is  conserved.  This  conservative  scattering  approach  can 
lead  to  an  overestimate  of  radiance. 

In  order  to  provide  a  more  realistic  simulation  of  radiation  in  spectral 
regions  and  along  atmospheric  paths  where  multiple  scattering  (MS)  is  a  sig¬ 
nificant  contribution  to  the  source  function,  an  efficient  and  accurate  scat¬ 
tering  parameterization  has  been  incorporated  into  the  LOWTRAN  and  FASCODE 
models  • 
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A  two-year  research  program  supported  by  AFGL/OPI  was  undertaken  at 
Atmospheric  and  Environmental  Research,  Inc.  (AER).  The  overall  program 
objective  was  to  incorporate  a  uniform,  compatible,  accurate  and  efficient 
multiple  scattering  scheme  into  the  present  LOWTRAN/FASCODE  models  so  that  the 
multiple  scattering  effects  due  to  particulates  could  be  realistically 
treated.  The  first  phase  of  the  program  was  to  study  and  evaluate  the 
available  multiple-scattering  parameterizations  and  recommend  the  most  appro¬ 
priate  one  for  its  incorporation  into  the  transmittance/radiance  models. 
Constraints  of  efficiency,  accuracy,  and  code  structure  of  LOWTRAN/FASCODE 
were  the  key  considerations  of  the  evaluation  processes.  The  results  of  a 
trade-off  analysis  to  determine  an  appropriate  approach  are  summarized  in 
Appendix  A.  The  second  phase  of  the  study  was  to  implement  the  most  suitable 
multiple  scattering  parameterization  into  the  LOWTRAN  and  FASCODE  models,  test 
the  resulting  multiple  scattering  versions  and  document  the  codes. 

This  report  summarizes  the  work  accomplished  in  the  second  phase  of  the 
overall  effort.  Following  this  introduction,  the  selected  approach  for 
multiple  scattering  is  described  in  detail  (Section  2).  In  Section  3,  the 
implementation  of  the  multiple  scattering  parameterization  within  the 
LOWTRAN/FASCODE  models  is  discussed  in  detail.  Comparisons  of  the  implemented 
code  are  made  to  exact  treatments  and  to  one  another  in  Section  4.  Finally, 
recommendations  are  provided  in  Section  5  for  enhancing  the  multiple 
scattering  capability  provided  here. 
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2.  MULTIPLE  SCATTERING  APPROACH 


Selection  of  an  appropriate  treatment  of  multiple  scattering  (MS)  for 
application  to  the  LOWTRAN  and  FASCODE  models  is  severely  constrained  by 
competing  requirements  of  desired  efficiency  and  accuracy,  and  limitations 
imposed  by  the  inherent  code  structures  of  these  models.  Additionally,  it  is 
desired  to  provide  an  approach  which  is  uniformly  applicable  to  all  spectral 
regions  considered  and  equally  appropriate  for  implementation  within  both 
LOWTRAN  and  FASCODE.  Based  on  these  considerations,  the  MS  parameterization 
selected  consisted  of  a  finite  stream  approach  (using  two  streams  for 
simplicity)  to  approximate  the  scattering  source  function. 


This  approach  could  be  implemented  directly  in  FASCODE  since  it  is  essen¬ 
tially  a  monochromatic  calculation.  It  is  known,  however,  that  there  exist 
difficulties  in  calculating  the  transmittance/radiance  averaged  over  a  finite 
spectral  interval  in  a  nongray  gaseous  absorber  with  multiple  scattering 
because  the  commonly-used  band  models  are  not  applicable  (cf.  Stephens, 

1984).  Such  is  the  case  for  implementing  an  MS  treatment  within  LOWTRAN.  The 
best  approach  to  solve  this  problem  is  the  use  of  the  k-distribution  method, 
which  decouples  the  multiple  scattering  from  the  gaseous  spectral  integration 
so  that  the  available  (monochromatic)  multiple  scattering  algorithms  can  be 
used  directly.  For  LOWTRAN,  the  stream  approximation  is  performed  through  an 
interface  routine  consisting  of  the  k-distribution  method.  For  practical  pur¬ 
poses,  this  consists  of  decomposing  the  band  model  determined  optical  proper¬ 
ties  into  a  set  of  equivalent  monochromatic  calculations  which  are  then  summed 
to  give  the  spectrally  averaged  results.  The  configuration  of  the  multiple 
scattering  approach  for  LOWTRAN  and  FASCODE  is  summarized  in  Figure  2-1. 


As  mentioned  earlier,  the  multiple  scattering  parameterization  has  to  be 
accommodated  by  the  present  LOWTRAN  and  FASCODE  code  structures.  This  con¬ 
straint  is  particularly  important  for  FASCODE.  In  the  FASCODE  application, 
gaseous  absorption  is  evaluated  directly  from  the  line-by-line  calculation. 
Fluxes  required  for  the  stream  approximation  are  calculated  via  the 
parameterized  adding  method.  The  adding  method  is  particularly  consistent 
with  the  code  structure  of  both  radiance/transmittance  models  since  they  treat 
one  layer  at  a  time. 
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Fig.  2-1.  AFGL  transmittance/radiance  code  structure  —  present  models  (solid 
line)  with  MS  parameterization  (dashed  line). 


In  the  FASCODE  model,  for  example,  the  evaluation  of  layer  optical 
properties  always  commences  with  that  level  in  the  selected  path  with  the 
highest  pressure  and  the  selected  spectral  sampling  interval  decreases  with 
pressure.  This  approach  insures  that  the  layer  spectral  resolution  is 
consistent  with  the  decrease  of  Voight  line  widths  at  higher  altitudes.  From 
the  perspective  of  the  line-by-line  calculation,  this  method  is  computation¬ 
ally  quite  efficient.  However,  it  is  inconsistent  with  monochromatic  multiple 
scattering  treatments.  In  order  to  accomodate  the  sampling  requirements  of 
both  the  line-by-line  and  multiple  scattering  calculations,  the  FASCODE 
implementation  employs  the  adding  formalism  to  aid  in  the  merging  of  the 
scattered  fluxeb  from  successive  layers. 

2.2  Theory 

We  provide  a  generic  outline  of  the  basic  theory  here.  This  prescription 
is  modified  slightly  for  specific  application  to  the  LOWTRAN  and  FASCODE 
models.  Details  of  the  specific  model  implementations  are  described  in 
Section  3. 


2.2.1  Radiance  and  Source  Function 

The  desired  radiance  1  at  wavenumber  v  for  an  arbitrary  path  with 
zenith  angle  cosine  and  azimuth  angle  (u,$)  is  given  by  the  solution  to  the 
radiative  transfer  equation  (RTE): 


(2.  1) 


U  Ly  (T  ,U  ,4>  )  =  Ly  (t  ,U  ,4>  )  -  Jv(t  ,\i  ,$  ). 

Here  is  the  optical  thickness,  y  is  then  cosine  of  the  path  zenith  angle, 
end  $  is  the  azimuth  angle  relative  to  the  sun’s  azimuth.  Vertical  optical 
depth,  ty,  will  depend  on  the  relevant  mechanisms  determining  the  extinction 
of  electromagnetic  radiation  for  the  spectral  region  characterized  by 
wavenumber,  v.  In  general  these  mechanisms  include:  (a)  molecular 
absorption,  ka,  (b)  molecular  scattering,  kg,  (c)  particulate  absorption,  oa. 
and  (d)  particulate  scattering,  ag.  Optical  depth  is  given  by  integrating  the 
relevant  vertical  extinction  coefficient  profiles  according  to: 

oo 

t  (z)  =  /  (k  (z)  +  k  (z)  +  o  (z)  +  o  (z)]  dz 
o  J  a  s  a  s 

z 

The  general  source  function,  Jy ,  including  scattering  of  solar  radiation  and 
thermal  emission,  is  given  by: 

J(r,y,<(>)  n  d0(T,y,<j>)  +  Jms(t,u,4i)  (2.2) 


where 


J0(T  »♦  ) 


-  nFe^o  p(n;-a0) 

+  [  l-u0(T  )]  B [G  (x  )  ] 


(2.3) 


and 


JmsCt  ,y  ,♦)  P(0;n'>  Kx.G'jdn'.  (2.4) 

Here,  wQ  is  the  single  scattering  albedo,  P  is  the  appropriate  angular 
scattering  or  phase  function,  and  B  is  the  Planck  function  at  temperature  0. 
The  extraterrestrial  solar  irradiance  is  given  by  F  and  the  path  and  solar 
directions  are  given  by  0  and  fl0,  respectively.  The  first  terra  in  Equation 
(2,3)  is  the  single  scattering  of  solar  radiation  while  the  second  is  the 
local  thermal  emission. 
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Radiance  solutions  to  the  R'E  (Eq.  (2.1))  are  subject  to  boundary  condi¬ 
tions  at  the  top  of  the  atmosphere  (t  *  0.0)  for  downward  radiance  and  at  the 
earth's  surface  (t  =*  x*)  for  upward  radiance.  At  x  *  0.0,  downward  diffuse 
radiance  from  space  is  zero  (the  direct  solar  irradiance  is  accounted  for  via 
the  primary  source  function)  resulting  in: 


1^(0 ,-y  ,4>  >  =  0.0 


(2.5) 


(At  millimeter  wave  frequencies  a  contribution  due  to  emission  at  the  cosmic 
background  effective  temperature  of  2.7  K  may  be  included  when  high  accuracy 
is  required.) 

Boundary  conditions  at  the  surface  will  depend  on  the  nature  of  the  sur¬ 
face  ref lectance/emittance  properties.  The  most  common  assumption  is  that  of 
Lambert  reflectance,  i.e.  the  upward  isotropic  flux  given  by  a  constant  sur¬ 
face  albedo,  r,  times  the  downward  flux.  Upward  and  downward  fluxes  F*( x)  at 
optical  depth  t,  are  defined  respectively  as: 

2n  1 

F  (t  )  «  /  f  1(t  ,±y  ,$ )ydyd(j>  (2.6) 

o  0 

This  results  in  a  lower  boundary  condition  upward  radiance  of: 

2n  1 

lb(i*,jj  ,4> )  *  -  [irFpoexp(-T*/uo)  +/  /  I(x*,-y ,$)ydyd$]  +  (1-r)  BlTs ]  (2.7) 

o  o 


The  three  terms  on  the  R.H.S.  of  (2.7)  are  respectively:  (a)  reflectance  of 
attenuated  solar  irradiance  (in  UV,  visible,  near  IR  spectral  regions),  (b) 
reflectance  of  downward  scattered  radiance  field,  and  (c)  thermal  emission  due 
to  the  surface  at  temperature,  Ts.  The  surface  emissivity  is  unity  minus  the 
surface  albedo,  i.e.,  (1-r).  LOWTRAN  requires  that  the  surface  albedo  be 
specified,  while  FASCODE  asks  for  the  surface  emissivity. 

General  radiance  solutions  to  the  RTE  for  upward  and  downward  radiances, 
respectively,  are: 


1  (t  ,+u  ,<p ) 
I  (t  ,-u  ,<j>  ) 


=  lb(T*,U,<fr) 

“  I fa( o ,  — H  ,4» ) 


-(t*-t)/u  .  s  -(t-x)/u  dt 

»  +  j  J(t,u,4>)  e  — 

i  ^ 

-x/y.  fT  .x  -(x-t)/y  dt 

i  +  J  J(t,y,$)  e  — 


(2.8) 

(2.9) 
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where  Che  1^  are  given  by  Che  boundary  condiCions  (2.7)  and  (2.5)  above, 
respecCively .  Incorporacing  chese  boundary  condiCions,  che  radiance  soludons 
become : 


1  2jt 


A  /  A  *•  II 

I(t,p,4»)  ■  [iTuoFe-T  yo  +  /  /  Kx*,— u*d>) wdud4> ] 

0  °T* 

+  (1-r)  B [T(t * )}  exp-(r*-T)/u  +  /  J(c,p,<$>)  e  T^M  —  (2.10) 


I  (t  ,  — u  ,4> )  =  J  J(c,u,4>)  e_(T_t)/u  — 


(2.11) 


In  Che  scream  approximaCion,  che  mulciple  scactering  conCribuCion  Co  Che 
source  funcCion  (Eq.  (2.4)  above)  is  approximaced  by  assuming  consCanc 
scatcered  radiances  I+  and  I  over  upward  (0+)  and  downward  (fi  )  hemispheres, 
respectively,  or  from  Eq.  (2.4): 


U+<t)  /  P(n,n+)da* 


+  J  p(fl,n“)dn"]. 


(2.12) 


Integrating  over  the  angular  scattering  functions  for  the  resulting  azimuthal- 
ly  averaged  backscatter  fractions,  0(u),  as  a  function  of  zenith  angle  cosine 
and  substituting  the  corresponding  fluxes: 


I*(t )  *  F±(t )/u 


(2.13) 


results  in 


jMS^t  »*» A )  = 


o)0(t  ) 


{f±(t)[1-B(u>I  +  F+(t)  0(u) }  < 


(2.14) 


This  simple  expression  for  the  multiple  scattered  contribution  to  the  source 
function  is  added  to  the  single  scattering  and  thermal  emission  contributions 
for  the  general  source  function,  J  (Eq.  (2.2)).  The  source  function  is  then 
integrated  along  with  the  desired  path  (as  in  Eqs.  (2.10)  or  (2.11))  to  obtain 
the  desired  total  radiance  including  the  approximated  MS  contribution. 
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The  evaluation  of  the  approximated  MS  source  function  (Eq.  (2.14)) 
requires  local  fluxes  F+,  F~,  backscatter  fractions,  B(u),  and  single 
scattering  albedos,  co0.  The  backscatter  fractions  ate  given  as  functions  of 
zenith  angle  cosine  and  asymmetry  factor  by  Wiscombe  and  Grams  (1976)  (see 
Figure  2-2).  A  small  error  is  introduced  by  assuming  these  backscatter 
fractions  for  the  equivalent  Henyey-Greenstein  phase  function  rather  than 
integrating  the  actual  function.  The  single  scattering  albedo,  to0(T),  for  a 
given  layer  with  total  optical  thickness  At  is: 


<d0(t)  -  At8/At 


(2.15) 


where  Aig  is  the  total  scattering  optical  thickness  of  the  layer.  Discreti¬ 
zing  equation  (2.14)  for  a  given  layer,  N,  the  contribution  of  multiple 
scattering  is  approximated  as: 


w0N 


JsX  “  Jo  +  {FN±[l-0N(u,gN)] 

+  FN?SN(u,gN)}. 


(2.16) 


Here  the  fluxes  are  taken  as  the  layer  mean  quantities  evaluated  at  a  level 
halfway  through  the  layer.  The  asymmetry  factor,  g,  is  a  measure  of  the 
directional  scattering  and  can  be  evaluated  from  the  phase  function. 

Once  the  source  function  is  approximated,  the  path  radiance  can  be 
evaluated.  Along  a  path  consisting  of  layers  (N)  and  the  layer  above  (N+l) 
with  transmissions  Tjj  and  Tjg+^,  respectively,  for  example,  the  emission,  E, 
depends  on  the  path  integral  of  the  total  source  function: 


Ci  -  Cn*i  +  ■'sANtl<+“X1-IN+i) 


(2.17) 


for  downward  looking  and 


Ci  ■  f,n  +  0-W 


(2.18) 


for  upward  looking,  where  the  intrinsic  layer  emission  is 


(1  "V  JSA(±^)- 


(2.19) 
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Figure  2-2.  Backscattered  fractions  5  and  B(y)  for  the  Henyey-Greenstein 
phase  function  versus  the  asymmetry  factor  g  for  a  range  of 
values  of  y  (Wiscorabe  and  Grams,  1976). 


2.2.2  Layer  Fluxes 
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Fluxes  approximate  the  required  radiances  for  evaluation  of  the 
multiply  scattered  source  function.  Upward  and  downward  fluxes  (F+  and  F~, 
respectively)  for  Individual  isolated  layers  are  evaluated  using  an 
appropriate  flux  parameterization.  For  example,  for  solar  scattering,  the 
hybrid  modified  delta  eddington  approximation  (Meador  and  Weaver,  1980)  is 
used.  The  chosen  flux  parameterization  also  provides  intrinsic  layer 
reflection  and  transmission  functions,  R  and  T.  These  fluxes  are  calculated 
using  standard  two  stream  parameterization  approaches.  To  accommodate  the 
flux  parameterizations ,  optical  properties  for  the  whole  atmosphere  (i.e., 
surface  to  space)  are  required.  The  approach  for  calculating  fluxes  thus 
consists  of  two  steps:  (1)  evaluating  local  layer  (i.e.,  intrinsic)  fluxes 
for  each  atmospheric  layer,  and  (2)  combining  these  to  obtain  the  actual  flux 
profiles  using  the  adding  method. 

Upward  and  downward  layer  fluxes  for  solar  radiation  are  given  by: 


_+  .  kT  .  _  -kx  .  _  -t/u 

F  ■  Ae  +  Be  +  Ce  o 


F“  {A(Yl-k)ekT  +  B(Y1+k)e"kT  +  Ye~T/uo} 


where  the  appropriate  constants  are  given  by: 

2  2 
k  -  (Yf  “  Y2  ^ 


(2.20) 

(2.21) 
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A  -  [B(Yx+k)  +  Y ) / C k-y x > 


„  kx*  .  „  -x */p  w,_  kx*  .  _  -kx*.. 

B  ■  (E^e  +  o)/(E^e  +  E^e  ) 

B(m  )  W* 

c  "  ™0t- - Tib(uo>  "  Y2U-e(yo)  1K-— — j) 

o  1-k  y 

o 

Y  *  C(y  ,  +  — )  ”  xFu  B(u  ) 

1  y  o  o 
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(2.22) 
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Ej  -  Y[l/(Yl-k)  -  r/Y2l 


[  -C  +  it  Fy  r  +  — 1 


E3  =  (Yj  +  lOU/dt-Yj)  +  r/Y2] 


[1  -  r(y j  +  k)/Y2l  (2.23) 

-{ l-g2-uo(4-3g)-wog2(4Bo+3g-4)} 

4[l-g2(l-U0)] 

{ 7-3g2-uo<4+3g)+uog2(4Bo+3g) } 

4[l-g2<l-M0)] 


Here,  r  is  the  Lambert  surface  albedo  and  the  solar  zenith  angle  cosine  is  uc 

The  transmission  and  reflection  functions  used  later  in  the  flux  adding  are 
given  by: 


R  -  F+/u0itF 
T  «  F~/u  ttF  +  exp 


( — r  /u  ) 


(2.24) 


For  the  thermal  fluxes,  a  linear  Planck  function  relation  across  an  atmo¬ 
spheric  layer  is  used.  In  so  doing,  the  parameterized  two-stream  solutions 
for  emission  from  the  layer  top  and  layer  bottom,  and  for  total  transmission 
and  reflection  are 

F+  »  a(PBt  -  mQ  -  Bfa)/D 

F"  -  a(PBb  +  mQ  -  Bt)/D 

T  »  a/D  (2.23) 

R  -  uv(eT 1  -  eT 1 ) /D 
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where  Bt  and  Bb  are  the  Planck  intensity  at  the  layer  top  and  bottom  and 
a2  -  1  -  uj0 
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(Bfe  -  Bt)/t 
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m 

ve  1  +  ue”11 
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ve1 1  -  ue-1 1  -  a 

(2.26) 
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The  optical  thickness  t  and  the  single  scattering  albedo  (dq  are  given  by 


t  -  ku  +  Tg(l  -  g)  +  Ta 
«0  "  Ts(i  "  g)/T 


(2.27) 


Here  k  is  the  gas  absorption  coefficient  (for  a  particular  wavelength  and 
probability  interval),  u  is  the  gas  amount,  x„  is  the  scattering  optical 
thickness,  x_  is  the  absorption  optical  thickness  and  g  is  the  asymmetry 
factor  for  the  particulate  matter  in  the  layer. 


2.2.3  Flux  Adding  Method 

To  obtain  the  actual  flux  profile  throughout  the  atmosphere,  intrinsic 
layer  fluxes  are  combined  algebraically  using  the  adding  method.  In  this 
method,  fluxes,  reflections,  and  transmissions  are  used  to  add  individual 
layers  together.  Composite  upward  fluxes,  and  reflection  functions  Rjjj, 

obtained  upon  adding  two  isolated  layers,  N  and  (N-l)  are  given  by: 


V-  FN+Tn(  1  FN-1+FNRN-1 ) ( 1-RNRN-1  r1  <2‘28> 

rN  “  RN  +  RN-1T&  (l  "  RNRN-i)  1 •  (2.29) 
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Analogous  equations  provide  composite  downward  fluxes  and  reflection  functions, 
and  R^,  respectively.  The  composite  upward  and  downward  fluxes  provide 
the  actual  upward  and  downward  fluxes  at  layer  interfaces  including  the 
effects  of  all  layers  above  and  below.  For  example,  the  upward  and  downward 
fluxes  at  the  boundary  between  layers  N  and  (N+l)  are  given  by: 

2F;  -  (iF;  +  ^rJhi  -  rJ%+i p.  (2.30) 

2pN+l  *  (lpN+l+1  FNRN+1^  ^^n+i)  (2.31) 

Once  obtained,  these  fluxes  are  substituted  into  Eq.  (2,16)  above  to  provide 
the  approximation  of  the  MS  source  function. 

2.2.4  Band  model  considerations 

For  LOWTRAN,  it  is  necessary  to  integrate  the  spectral  radiance  values 
over  a  finite  spectral  interval  (~  20  cm-1).  The  basic  problem  encountered  in 
the  calculation  of  radiative  transfer  in  low  spectral  resolution  in  a  hazy  or 
cloudy  atmosphere  is  the  coupling  between  the  processes  of  scattering  and  ab¬ 
sorption  due  to  cloud/aerosol  particles  and  absorption  by  atmospheric  gases. 
The  main  difficulty  is  that  the  integration  over  frequency  cannot  be  properly 
accounted  for  by  the  usual  band  model  technique  for  gaseous  absorption  because 
they  do  not  allow  for  multiple-scattering.  A  direct  line-by-line  integration 
over  frequency  would  be  very  time  consuming.  One  alternative  way  of  carrying 
out  the  frequency  integration  is  to  use  the  "k-distribution  method"  for  homo¬ 
geneous  layers  (Arking  and  Grossman,  1972)  and  the  "correlated  k-distribution 
approximation"  for  inhomogeneous  atmospheres  (cf.  Wang  and  Ryan,  1983). 

For  gaseous  absorption,  the  k-distribution  method  is  comparable  to  line- 
by-line  calculations  (Arking  and  Grossman,  1972),  This  method  is  equivalent 
to  the  exponential-sum  fitting  method  (see  Wiscombe  and  Evans,  1977)  and  to 
the  path  length  distribution  method  (see  Bakan  et  al.,  1978).  However,  in 
general,  the  latter  two  methods  use  scaling  approximations  to  account  for 
atmospheric  inhomogeneity  while  the  correlated-k  approximation  assumes  certain 
relationship  between  k  values  at  different  pressure  and  temperature  levels. 

The  accuracy  of  the  approximation  is  excellent  for  the  9.6  pm  0-j  band  thermal 
radiation  calculations  (see  Lacis  et  al.,  1979). 


Yamamoto  et  al.  (1970,  1971)  used  finite  sums  of  exponentials  to  describe 
the  non-gray  nature  of  water  vapor  absorption  and  carried  out  solutions  of  the 
equation  of  transfer  for  homogeneous  band  layers  using  both  Chandrasekhar's 
principles  of  invariance  as  well  as  the  discrete  ordinate  technique.  Both 
techniques  require  extensive  numerical  calculations.  On  the  other  hand,  two- 
stream  approximations  together  with  the  correlated-k  approximation  have  been 
used  to  study  the  radiative  effects  of  aerosols  (see  Hansen  et  al.,  1980).  As 
best  summarized  recently  by  Stephens  on  efficient  and  accurate  radiation 
parameterizations  (1984,  p.  862  of  the  paper):  "...  only  the  k-distribution 
approach  can  be  readily  incorporated  into  scattering  models. . . " 

In  a  homogenous  gas  layer,  the  k-distribution  function  is  formerly 
related  to  the  mean  transmission  function  T^v(u), 


TavCu)  2  ^  I  e‘k“  dv  5  /  £‘k>*"k“ 


Av 


dk 


Av 


. 1  -ku  n  ~kiu 

=  J  e  dg  **  I  e  Ag^ 
i*l 


(2.32) 


where  Av  is  the  narrow  repeated  interval  (20  cm"1  in  LOWTRAN)  and  u  is  the  gas 
amount.  The  f(k)  for  a  given  gas  at  a  specified  Av  is  the  probability  density 
function  such  that  f(k)dk  is  the  fraction  of  the  frequency  interval  for  which 
the  absorption  coefficient  is  between  k  and  k+dk.  Eq.  (2.32)  reveals  that  the 
transmission  depends  on  the  distribution  of  k-values  within  Av,  but  not  on  the 
ordering  of  the  values.  The  cumulative  k-distribution  function  is  g(k),  while 
( k ^ ,  Ag^)  are  the  discreet  sets  of  values  to  approximate  the  integral. 

By  expressing  the  band  model  transmission  as  the  sum  of  exponentials,  the 
multiple  scattering  calculation  for  each  component  can  be  performed 
independently  as  if  it  were  a  monochromatic  problem.  These  are  weighted  and 
summed  (as  in  Equation  (2.32)  to  recover  the  essential  band  model  character  of 
the  problem. 
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The  fit  of  Wlscombe  and  Evans  (1977)  has  been  used  for  the  two  LOWTRAN 
transmission  functions  of  water  vapor/unlformly  mixed  gases,  and  ozone.  The 
accuracy  of  the  fitting  Is  In  general  within  a  few  percent  for  T  >  0.1.  For 
Inhomogeneous  atmospheres,  we  adopt  the  same  scaling  approximation  used  in 
LOWTRAN,  l.e. , 


ki(P,0) 


kf(  Pq  ^ 


(2.33) 


where  0O,  PQ  are  reference  temperatures  and  pressures,  respectively. 


-19- 


3.  IMPLEMENTATION  OF  MULTIPLE  SCATTERING  PARAMETERIZATION 


3.1  FASCODE  Implementation  of  the  Stream  Approximation  for  Multiple  Scattering 

3.1.1  Overview 

The  existing  code  structure  of  the  FASCODE  line-by-line  algorithm  pre¬ 
sents  a  constraint  for  the  implementation  of  a  multiple  scattering  capability. 
In  order  to  optimize  computational  efficiency  and  minimize  online  storage  re¬ 
quirements,  FASCODE  program  logic  calculates  path  transmittance  and  radiance 
In  a  single  unidirectional  pass  through  the  desired  atmospheric  path.  This 
pass  commences  with  the  highest  pressure  sublayer  (i.e.,  that  requiring  the 
least  spectral  resolution)  and  progresses  monotonically  to  the  lowest  pressure 
sublayer  requiring  the  highest  spectral  resolution.  Cumulative  quantities 
with  spectral  structure  (such  as  transmittance  and  radiance)  are  evaluated  by 
merging  the  lower  resolution  results  from  the  previous  layer  with  the  higher 
resolution  results  of  the  current  layer.  While  this  procedure  is  appropriate 
for  nonscattered  radiance  where  the  local  source  function  for  a  layer  consists 
solely  of  contributions  from  local  thermal  emission,  the  methodology  is  not 
valid  in  the  presence  of  multiple  scattering.  Fundamentally,  the  local  mul¬ 
tiply-scattered  (MS)  source  function  is  composed  of  contributions  scattered 
from  layers,  both  above  and  below  the  current  layer  In  addition  to  local 
thermal  emission.  Therefore,  a  modification  of  the  FASCODE  program  logic  is 
required  to  Incorporate  the  stream  approximation  (SA)  for  multiple  scattering 
within  FASCODE.  The  SA  efficiently  provides  an  estimate  of  the  local 
multiply-scattered  source  function  given  the  availability  of  a  suitable 
estimate  of  the  profiles  of  upward  and  downward  flux.  Evaluation  of  the  flux 
profile,  however,  requires  a  priori  total  atmospheric  properties  such  as 
optical  depth  and  single  scattering  albedo  profiles.  This  is  not  possible 
with  the  current  layer  by  layer  FASCODE  algorithm.  Once  obtained,  the  source 
function  for  the  multiply-scattered  radiance  may  be  summed  along  the  desired 
path  in  analogy  to  the  approach  currently  used  within  FASCODE  to  obtain  the 
thermal  emission  via  summation  of  the  local  Planck  function  contributions. 

This  apparent  inconsistency  of  FASCODE  with  multiple  scattering  is 
resolved  by  combining  the  stream  approximation  approach  with  the  adding  method 
to  obtain  the  required  MS  flux  profiles.  In  the  adding  method,  properties  of 
individual  layers  (such  as  emission,  transmission,  and  reflection)  are 
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combined  algebraically  to  yield  the  overall  properties  of  the  combined 
layers.  A  typical  application  of  the  adding  method  consists  of  two  opposite 
unidirectional  passes  throughout  the  atmosphere.  In  the  FASCODE 
implementation,  the  existing  upward  pass  is  exploited  to  set  up  one  part  of 
the  adding  calculation.  This  is  done  subsequent  to  the  existing  calculation 
of  optical  depths  via  the  combined  line-by-line  gas  absorption  and  LOWTRAN- 
based  aerosol  extinction  routines.  The  upward  pass  also  provides  an 
opportunity  to  evaluate  and  merge  radiance  contributions  due  to  thermal 
emission  only.  A  downward  pass,  unique  to  MS  cases,  completes  the 
calculation.  During  the  downward  pass,  actual  fluxes  will  be  calculated  for 
each  layer,  layer  multiply-scattered  source  function  evaluated  using  the 
stream  approximation,  and  multiply-scattered  radiances  summed  along  the  de¬ 
sired  path.  These  radiances  will  be  merged  with  the  nonscattered  radiance 
contributions  evaluated  during  the  upward  pass  to  yield  total  cumulative  mul¬ 
tiply-scattered  radiance. 

The  combined  stream  approximation/adding  method  approach  for  FASCODE  is 
summarized  in  the  flow  chart  presented  in  Figure  3—1.  The  shaded  areas  denote 
modifications  to  FASCODE  to  accommodate  multiple  scattering.  The  unshaded 
code  segment  on  the  left-hand  aide  of  the  diagram  represents  the  existing  non- 
MS  FASCODE  algorithm.  Within  the  shaded  portion  during  the  upward  pass  are 
code  modifications  to  calculate:  (a)  layer  single  scatter  albedo,  w0N,  and 
phase  function  asymmetry  factor,  g^;  (b)  layer  reflection  and  transmission 

functions,  RN,  TN,  respectively  and  intrinsic  upward  and  downward  fluxes  FN* 
using  the  two  stream  approximation;  (c)  beginning  of  the  adding  method  for 
composite  upward  fluxes  and  reflection,  *F+,  R+,  respectively;  and  (d)  layer 
and  merged  values  of  the  thermal  radiance  contribution  EN°  based  on  the  JN° 
source  function  described  above.  The  downward  pass  completes  the  adding 
method.  Details  will  be  described  in  the  following  sections.  Notably,  the 
time-consuming  line-by-line  calculations  are  done  once  only  during  the  upward 
pass.  The  downward  pass  is  a  simple  merging  of  precalculated  and  stored 
quantities  based  on  algebraic  manipulations. 
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3.1,2  Adding  Method 


To  illustrate  the  application  of  the  adding  method  to  calculate  the 
actual  flux  profiles,  consider  the  layered  atmosphere  illustrated  in  Figure 


3-2.  There  are  four  layers  ( N— 1 ) ,  (N),  (N+l),  and  (N+2)  above  an  emitting 

surface  (N-2).  As  each  layer  is  treated  during  the  existing  FASCODE  upward 

pass,  the  following  quantities  are  evaluated: 

N+2 

TN+2*  TN+2’  FN+2 

N+l 

RN  + 1 ’  TN+ 1  *  FN+1 

N 

rn*  tn*  fn 

N-l 

rn-i*  tn-i*  fn-i 

N-2 

///////////////  Rn_2.  F^j-2 

Figure  3-2.  Schematic  of  FASCODE  model  atmosphere  layers. 
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gas  absorption  optical  depth  via  line-by-line 

molecular  scattering  via  LOWFAS 

aerosol  absorption  optical  depth  via  LOWFAS 

aerosol  scattering  optical  depth  via  LOWFAS 

mean  layer  thermodynamic  temperature  via  FASCODE  layering. 


For  a  multiple  scattering  (MS)  calculation,  the  shaded  branch  of  the  upward 
path  in  Figure  3-1  is  selected  and  the  following  two  additional  quantities  are 
calculated: 

v  N 

w0‘  single  scattering  albedo  =  (t^  g 
g'  phase  function  asymmetry  factor 

The  asymmetry  factor  is  obtained  from  a  look.-up  table  for  the  appropriate 
aerosol  model,  wavenumber  domain,  and  level. 
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The  asymmetry  factor  is  obtained  from  a  look-up  table  for  the  appropriate 
aerosol  model,  wavenumber  domain,  and  level. 

Next,  a  simple  two-stream  approximation  is  used  to  evaluate  the  intrinsic 
layer  fluxes,  transmission  function  (direct  and  diffuse),  TN,  and 

reflection  function,  R^.  The  intrinsic  layer  fluxes,  FN±  are  the  upward  flux 
at  the  top  and  the  downward  flux  from  the  bottom  of  the  Nth  isolated  layer, 
respectively.  These  are  analytical  transformations  of  the  form: 


,  N  N  N  N  N, 
IT  .  t  ,  x  ,  g  ,  8 
1  g  a, a  a, s’  6  J 


iFN  *  FN  ’  V  V' 


(3.1) 


Specific  two  stream  flux  formulations  for  solar  and  thermal  scattering  are 
given  in  Section  2,2.2.  These  are  evaluated  once  during  the  upward  pass 
through  FASCODE  at  increasing  spectral  resolution  with  decreasing  layer  mean 
pressure,  i.e.,  Av(N+l)  <  Av  (N)  where  N  is  layer  index  numbered  consecutively 
from  surface  to  space,  and  Av  is  the  spectral  resolution  for  each  layer 
determined  by  the  FASCODE  algorithm. 

Surface  emission  is  treated  as  the  zeroth  layer  with: 


Fo+  -  iresB(0s) 

Ro  “  U  -  es)  (3.2) 


where  8g  is  the  surface  temperature  and  eg  is  the  surface  emissivity.  Like¬ 
wise,  a  target  at  the  top  of  the  path  may  be  treated  as  the  Mth  layer  with: 


I’m  “  ttc  tBC©  c  ) 
Rm  =  (1  “  et) 


(3.3) 


Similar  considerations  can  be  applied  to  treat  nonscattering  layers  above  and 
below  the  scattering  layer  system.  As  noted  in  the  flow  chart  (Figure  7  ), 

an  additional  modification  necessitated  in  order  to  accomplish  a  multipl 
scattering  calculacton  is  the  consideration  of  bounding  layers  above  ar  below 
the  altitudes  of  the  nodes  of  the  specified  path  (i.e.,  HI  and  H2).  A  minimum 
of  one  layer  extending  from  the  surface  to  HI  and  from  H2  to  13  km  (assuming 
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H2  <  15  km)  Is  Included.  In  order  to  Illustrate  application  of  the  adding 
method,  the  required  relations  for  the  system  illustrated  in  Figure  3-2  are 
presented  explicitly.  It  has  been  assumed  that  multiple  scattering  contribu¬ 
tions  above  15  km  will  be  negligible. 

3. 1.2.1  Upward  Adding  Pass 

The  loop  begins  by  adding  layer  (N-l)  to  layer  (N-2).  Subsequent 
layers  are  added  one  at  a  time  until  layer  (N+2)  is  reached.  The  modified 
upward  flux  from  the  combination  of  (N-2)  and  (N-l),  consists  of  three 

contributions  (see  Figure  3-2): 

i,  the  intrinsic  upward  flux  from  layer  (N-l),  FjV 

it.  the  multiple  reflection  of  the  intrinsic  upward  flux  from  layer 
(N-2),  Fjjj_2>  transmitted  through  layer  (N-l): 

1,e*’  FN-2  +  TN-1  +  FN-2  TN-1  Vl  V-2  + 

FN-2  Vl^Vl  ^-2^  +  *  *  * 

FN-2  Vl^1  “  RN-1  1 


ill,  the  multiple  reflection  of  the  intrinsic  downward  flux  from  layer 
(N-l),  F^_p  reflected  from  layer  (N-2)  and  transmitted  through 
layer  (N-l): 


i,e,»  fn-i  Vi  rn— 2  +  fn-i  tn-i  RN-2^RN-l  RN-2 ^  +  •  •  • 

FN-I  TN-1  ^-2^  “  RN-3.  V2) 

These  will  sum  to  the  composite  flux: 

lfN-l  =  FN-1  +  TN-lV-2  +  Vl  V2V  ”  Vl  V2^ 

Similar  considerations  lead  to  the  composite  reflection  of  layers  (N-l) 


(3.4) 


and 


(N-2): 
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(3.5) 


hi-l  "  *11-1  +  *N-2  TN-1^  “  *N-1 


-1 


Calculated  composite  layer  values  are  required  for  the  evaluation  of  the  com¬ 
posite  values  of  the  subsequent  layers.  These  intermediate  quantities  are 
merged  to  the  next  higher  layer's  spectral  resolution  and  stored  for 
subsequent  use  on  the  downward  pass.  At  the  termination  of  the  upward  pass, 
the  following  quantities  have  been  calculated  for  each  layer  and  stored  for 
uac  during  the  downward  pass:  (a)  the  Intrinsic  layer  values,  u>0N,  RN,  TN, 
(the  intrinsic  upward  fluxes  are  not  required),  (b)  the  composite  quantities, 
1  Fjg ,  R  ,  and  the  nonscattered  contribution  to  emission,  E^j.  These  values  are 
stored  only  at  local  layer  resolution  to  minimize  total  memory  requirements. 
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Next  layer  N  is  added  to  the  composite  of  (N-l,  N-2),  resulting  in 

composite  values  for  layer  N: 

y$ 

ife 

lpN  ”  FN  +  TN^FN-1+  PN  S-l^1  "  ^-1^ 

(3.6) 

& 

^  "  *N  +  hl-l  tn  ^ 

(3.7) 

V 

v  - 

:S 

V 

This  process  is  repeated  for  layer  (N+l): 

lpN+l  =  PN+1  +  TN+l(lpN  +  PN+1  RN^1  ”  ^+1  hi  ^ 

(3.8) 

% 

S+l  “  V+l  +  hi  TN4l^1  “  ^+1^  ) 

(3.9) 

k*i 

1 

*  • 

and  layer  (N+2): 

'{j 

1  *N+2  =  PN+2  +  TN+2^ 1 PN+1  +  FN+2  S+l^1  “  *N+2  *N+1  ^ 

(3.10) 

ri 

Iv 

S+2  “  ^+2  +  RN+l  TN+2^  "  *N+2  *14+1 ) 

(3.11) 

1 

Since  there  are  no  layers  above  (N+2),  the  upward  composite  *-8 

#-  ■ 

#  r 

identically  the  actual  upward  flux,  i.e., 

,v 

.y 

.  0 

2n  +  1_  + 

F  ss  p 

N+2  N+2 

(3.12) 
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3. 1.2. 2  Downward  Pass 


The  downward  pass  is  devoted  exclusively  to  the  adding  method 
and  stream  approximation.  No  additional  line-by-line  calculations  are 
performed.  In  analogy  to  its  upward  counterpart,  downward  composite  fluxes, 
and  reflectance  functions,  RN~,  are  calculated  during  the  unidirectional 
downward  pass.  However,  in  addition  both  actual  upward  and  downward  fluxes, 
2FjJ±,  are  calculated  from  the  composite  fluxes  and  reflection  functions. 

For  the  topmost  layer,  no  downward  composite  fluxes  are  possible  (unless 
there  are  nonscattering  layers  above  N+2).  Therefore,  the  actual  upward  and 
downward  fluxes  at  the  interface  of  layers  (N+l)  and  (N+2)  are: 


2FN+1  =  (lfN+l  +  FN+2  RN+1^1  “  *Sl+2  RN+l) 

2pN+2  *  (lfN+2  +  lfN+l  RN+2^1  ~  ^N+2  ^+1^ 


(3.13) 

(3.14) 


At  this  point,  actual  upward  and  downward  fluxes  are  available  at  the 
boundaries  of  layer  (N+2),  facilitating  calculation  of  the  multiply-scattered 
source  function,  3^2*  usin£>  the  stream  approximation  (see  Equation  (2.14)). 
This  source  function  provides  the  layer  contribution  to  multiply-scattered 
emission,  This  is  added  to  the  thermal  contribution  E^+2  evaluated 

during  the  upward  pass. 

Proceeding  downward,  the  composite  quantities  are  first  calculated  by 
adding  the  next  layer  down  to  the  composite  above  it,  for  example,  the 
downward  composite  flux  emerging  from  layer  (N+l)  and  the  corresponding 
downward  composite  reflection  are: 


lFN+l  ”  FN+1  +  TN+l^FN+2  +  FN+1  RN+2^1  ”  *N+2  S+l^ 
RN+1  =  RN+1  +  ^+2  W1  “  S+l  ^+2^ 


Actual  fluxes  at  the  interface  of  layers  N  and  (N+l)  are  then  obtained  from 
the  composite  values: 

2fJ  ■  +  ‘Vi  <)(>-<  v,r*  o-i 
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2iw  ■  I‘v, +  lFS  ’Will*  -  ^  V,r 


(3.18) 


For  che  next  layer,  N,  the  procedure  is  identical,  resulting  in  the  com¬ 
posites  : 


lFN  -  fn  +  tn^1fn+i  +  fn  *4^+1  H1  “  ^  ^+0 

^  ^  +  hi+i  V1  "  S  ^+1^ 


(3.19) 


(3.20) 


and  actual  fluxes: 


2pN-l  =  (1fn-i  +  lpN  “n-iH1  ”  ^-1  ^  ^ 
2FN  -  (lFN  +  lFN-i  h  Hi  -  Cl  h)~l 


(3.21) 


(3.22) 


Finally,  the  procedure  terminates  at  the  surface  layer  (N-2),  layer  ( N— 1 ) 
interfaces  where: 


-  v-1 


lpN-l  =  FN-1  +  TN-JlpN  +  FN-1  ^  H1  ^  ) 

RN-1  =  RN-1  +  RN  W1  ”  ^-1  S  ) 


(3.23) 


(3.24) 


2pN-2  ”  ^  FN-2  +  lpN-l  ^-2  Hi  "  ^-2  S-J 

2pN-l  "  (lFN-l  +  FN-2  VlH1  _  ^-2  RN-1  ^ 


(3.25) 


(3.26) 


In  this  last  step,  the  upward  surface  emission  flux,  F^_2>  and  reflection  func¬ 
tion,  R«4_2*  are  composite  quantities  since  there  are  no  layers  below.  This 
completes  the  downward  pass  and  the  multiply-scattered  radiance  calculation. 


3.1.3 


Lmatlon,  Source  iunct i on,  Radiance 


Using  the  stream  approximation  (SA),  the  multiply-scattered  source 
function  at  optical  depth,  t,  for  a  path  with  zenith  angle,  cos  !  u,  will  be 
given  by  (see  Equation  (2.16)): 
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Jc,(t,  ±g)  •  [1  -  a)  (t)]B<0)  +  -- —  {F±(t)[1  -  $(y)l  +  F  +  6(g)}  (3.27) 

The  plus  sign  is  for  downward  looking  paths  (u  >  0)  while  the  minus  sign  is 
for  upward  looking  paths  (p  <  0).  In  the  expression  above,  w0(t)  is  the 
single  scattering  albedo,  B(0)  is  the  Planck  function,  and  6(p)  is  the 
backscatter  fraction  for  path  zenith  angle  p.  The  backscatter  fraction  is  a 
function  of  the  particulate  asymmetry  factor  at  level  t,  g(x),  and  is 
available  from  Wiscorabe  and  Grams  (1976). 

For  level  N,  the  level  source  function  is  given  by: 

JNSA(±p)  -  d-(0oN)B  (eN)  + 

N 

”  {FN4[1  -  BN(p,gN)]  +  FN?BN(p,gN)}  (3.28) 


The  first  term  above  is  due  to  nonscattered  thermal  emission,  while  the  second 
terra  is  due  to  multiple  scattering.  Thus,  the  source  function  can  be  written  as 


J 


N 


SA 


(3.29) 


where  the  thermal  emission  terra  JN°  depends  only  on  local  layer  properties, 

MC 

and  the  multiply-scattered  term  JN  requires,  in  addition  to  local  proper¬ 
ties,  the  local  fluxes  which  depend  on  the  overall  path  properties.  The  first 
terra  is  a  generalization  of  the  nonscattering  thermal  emission  source  function 
currently  evaluated  within  FASCODE  and  may  be  treated  in  a  like  manner. 

The  path  radiance  or  emission  with  multiple-scattering  is  obtained  by 
summing  radiance  contributions  along  the  path  essentially  as  is  currently  done 
with  no  scattering.  Along  a  path  consisting  of  two  layers,  layer  N  and  the 
layer  above,  N  +  1,  the  emission  depends  on  the  path  integral  of  the  source 
function  and  is  given  by: 


JN+1 


en  tn+i 


N+l 


SA 


(+u)(l  -  w 


(3.30) 


for  downward  looking  and: 


(3.31) 


■'N+l 


+  (l  - 


TN+1J 


N+l 


SA 


<-p)T, 


for  upward  looking  where: 


U  -  TN)  4(±W. 


(3.32) 


With  no  scattering,  N+*  *  0  and  the  source  function  for  each  layer 

reduces  to  the  Planck  function  as  in  the  original  FASCODE,  i.e.,: 


lira  JoaN  =  8(e„).  (3.33) 

n  SA  N 

u0+  0 

As  described  above,  the  thermal  emission  source  can  be  summed  along  the 
path  in  analogy  to  the  current  nonscattering  radiance  calculation  to  provide 
the  path  thermal  emission  EN°.  Evaluation  of  the  multiply-scattered  radiance 
contribution  Ejj  for  each  layer  is  accomplished  via  the  adding  algorithm  on 
the  complementary  downward  pass  where  it  is  merged  with  the  previously  calcu¬ 
lated  nonscattered  contributions  to  obtain  total  radiance. 


3.1.4.  Description  of  FASC0DE2  with  Multiple  Scattering 

In  implementing  multiple  scattering  into  FASC0DE2,  the  main  considera¬ 
tions  were  in  maintaining  the  integrity  of  the  existing  program,  and  in  making 
the  resultant  code  as  transparent  to  the  user  as  possible.  The  first  and 
most  constraining  element  of  the  existing  code,  was  the  use  of  panels  in  doing 
the  radiance  calculation.  This  necessitated  that  any  multiple  scattering  im¬ 
plementation  also  be  adjusted  to  utilize  this  methodology. 
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All  of  the  routines  to  be  listed  below,  utilize  the  existing  FASCODE 
method  of  working  with  only  a  "panel"  of  numbers  at  a  time.  This  necessitates 
a  great  deal  of  1/0  for  doing  any  calculation  with  multiple  scattering,  but 
allows  for  complete  compatibility  with  the  existing  structure  and  execution  of 
FASCODE.  Another  benefit  of  this  approach,  is  that  the  program  is 
compartmentalized,  and  the  user  should  find  no  difference  in  run  time  for  this 
program  with  multiple  scattering  turned  off,  and  the  previous  version  of 
FASC0DE2 . 

As  discussed  above,  flux  adding,  used  in  conjunction  with  the  stream  ap¬ 
proximation  wa3  implemented  as  the  best  approach.  This  involved  the  extensive 


-30- 


.-  -a  —  j,  v 


v  r/  v  .v  w 


-V  V  V.V 


i  Y.Yi'-T.  Y 


modification  of  EMINIT,  and  RADMRG.  These  routines  were  modified  to  merge  the 
composite  fluxes  and  drive  the  flux-adding  during  the  radiance  calculation. 

The  primary  additions  to  the  code,  were  the  subroutines  FLXADD,  FLXDWN,  and 
SRCFCN.  FLXADD  calculates  and  adds  fluxes  for  each  of  the  FASCODE  layers. 

This  routine  is  called  within  the  layer  loop  from  both  EMINIT,  and  RADMRG. 
FLXDWN  drives  the  downward  flux  adding,  and  is  called  from  XLAYER  after  the 
FASCODE  upward  path  is  completed.  SCRFCN  evaluates  the  multiple  scattered 
source  function,  and  is  called  from  FLXDWN. 

Several  secondary  routines  that  were  added  were  MSEM1S,  TLEM1S,  XLAYMS, 
FLUXLP,  MSIN ,  and  MSOUT,  MSEMIS  merges  the  multiple  scattered  contribution 
to  the  radiance  for  each  layer.  TLEM1S  merges  the  total  multiple  scattered 
radiance  contribution  obtained  from  MSEMIS  and  merges  it  with  the  total  single 
scattered  radiance  calculated  within  FASCODE.  XLAYMS  drives  the  calculation 
of  the  radiance  contribution  from  layers  above  the  desired  path.  FLUXLP 
merges  the  fluxes  and  drives  the  flux-adding  for  layers  in  which  radiances  are 
not  evaluated.  The  purpose  of  these  last  two  routines  will  be  discussed  in 
Section  3.1.5.  MSIN  reads  merged  fluxes  from  the  previous  layer  while  MSOUT 
writes  out  the  merged  fluxes  for  the  current  layer. 

Some  minor  additions  were  the  functions  BETABS,  and  E3.  BETABS  evaluates 
the  backscattered  fraction  for  a  given  angle  and  asymmetry  factor.  E3  evalu¬ 
ates  the  exponential  integral  E3(x).  Some  minor  changes  were  also  made  to 
EMIN,  GETEXT,  AERF,  EXTDTA,  EXABIN ,  AEREXT,  TRANS,  LOWTRN,  and  ATMPTH .  EMIN 
was  modified  to  evaluate  the  thermal  emission,  for  multiple  scattering  cases. 
GETEXT  buffers  in  the  asymmetry  factors  along  with  the  extinction  and  scatter¬ 
ing  coefficients.  AERF  calculates  the  scattering  optical  thickness,  and  the 
asymmetry  factor  from  the  buffered  in  quantities.  EXTDTA  now  contains  the 
asymmetry  factors  for  each  of  the  LOWTRAN  models.  EXABIN  does  the  level 
weighting  for  the  asymmetry  factors.  AEREXT  evaluates  the  asymmetry  factors 
for  each  wavenumber.  TRANS  does  the  weighted  average  for  the  final  asymmetry 
factor  result.  LOWTRN  was  modified  to  read  in  the  asymmetry  factors  for  user- 
defined  aerosols.  ATMPTH  was  modified  to  determine  the  multiple  scattering 
path.  This  will  be  discussed  in  the  next  section. 

A  flowchart  for  FASC0DE2  with  multiple  scattering  is  illustrated  in  Table 
3-1.  The  structure  of  FASC0DE2  with  multiple  scattering  is  outlined  in  Table 
3-2  giving,  also,  brief  descriptions  of  the  altered  and  added  routines.  The 
segment  load  file  and  file  usage  are  shown  in  Tables  3-3  and  3-4,  respectively. 
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Table  3-1 

Flowchart  of  FAS CODE  Multiple  Scattering  Routines 

Underlined  subroutines  are  modified  from  the  existing  FASCODE; 
others  are  new  subroutines 
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Table  3-2 

Description  and  Segment  Load  Grouping  of  FASC0DE2  Subroutines 


(This  listing  will  contain  all  the  FASCODE  routines  but  will  only  describe 
the  routines  pertinent  to  the  multiple  scattering  calculation.) 


FASC0D2*  Fast  Atmospheric  Signature  Code 

/ADR1VE/ 

/CONSTS/ 

/FILHDR/ 

/IF1L  / 

/LANCHN/ 

/LAS IV  / 

/MAIN  / 

/MSACCT/*  Multiple  Scattering  files  and  parameters 

/MSCONS/*  Multiple  Scattering  flags  and  constants 

BUFIN 
BUFOUT 
COPYFL 
ENDFIL 
SK.IPFL 

FLTRFN 

CNVFLT 

RDSCAN 

LASER 

PLTFAS 

/AXISXY/ 

/NAME  / 

/PLTHDR/ 

/POINTS/ 

/T1TLOC/ 

/YCOM  / 

ENDPLT 

EXPT 

FSCLIN 

LINT 

MNMX 

TEMPFN 

TENLOG 

XNTLOG 

Legend:  /Common  Blocks/ 

(Block  Data) 

Subroutines 

NOTE:  designates  routines  added  or  modified  in  the  implementation  of 

Multiple  Scattering  by  AER. 
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Table  3-2  (Cont.) 


AX1SL 


AXLOG 


BBSCLE 


FPL1NE 


HEADER 


SCANFN 


CKPRNT 

CNVREC 

CONVSC 

INTER 

PANLSC 

RDSCAN 

SHAPEG 

SHAPET 

SHRKSC 

SINCSQ 


(BTEST) 


XLAYER* 


/ABSORA/ 
/ASYMMA/* 
/LINHDR/ 
/MLTSCT/* 
/SCATTA/ 
/RMRG  / 
/XBLF  / 
MS  COPY* 


Controls  Layer  by  Layer  Calculation 
Aerosol  Asymmetry  Factors 

General  Labeled  Common  for  Multiple  Scattering 


Copies  Multiple  Scattering  output  files  by  panels 


FLXDWN* 


SRCFCN* 


AERF* 

BBFN 

E3* 

FLXADD* 

GETEXT* 
MS  IN* 
MSOUT* 
RADFN 


BETABS* 

MSEMIS* 


Controls  Layer  by  Layer  Downward  Flux  Adding 
Evaluates  Extinction,  Scattering,  and  Asymmetry 
Computes  Black  Body  Function 
Computes  Exponential  Integral  E3 
Calculates  and  does  Adding  of  Fluxes  for  Current 
Layer 

Reads  Extinction,  Absorption,  and  Asymmetry 
Reads  Merged  Fluxes  from  Previous  Layer 
Writes  Merged  Fluxes  for  Current  Layer 
Computes  Radiation  Function 

Evaluates  Multiple  Scattering  Source  Function 
and  Radiance 

Calculates  Back-Scattering  Fraction 
Merges  Multiple  Scattered  Radiance 


,v.v 
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OPDPTH 

/ABSORB/ 
/SCATTR/ 
/LBLF  / 


CONTNM 

XINT 

FRNC02 

|  (BFC02) 

FRN296 

|  (BFH20) 

N2CONT 
|  (BN2) 

02C0NT 
|  (02C) 

03CHAP 
|  (03CH) 

O3HHT0 

(03HH0) 


03HHT1 

I  (03HHL) 


03HHT2 

I  (03HH2) 


SL296 

|  (BS296) 

SL260 

|  (BS260) 

HIRAC1 

I  /FNSH  / 

(BHIRAC) 
(BMOLEC) 
(BSHAPEL) 
( VOICON ) 
ABSOUT 
CNVFNV 
MOLEC 
PANEL 
QV 

RDLIN 

R1PRNT 

SHAPEG 

SHAPEL 

VERFN 


LBLF4 

(BMOLEC) 
(VO  ICON) 
CONVF4 
MOLEC 
QV 

RDLIN4 

SHRINK 

NONLTE 

/VBNLTE/ 

HIRACQ 

/ FNSQ  / 
(BHIRAQ) 
(BMOLEC) 
(BSHAPEL) 
(VOICON) 
ABSOUT 
CNVFNQ 
MOLEC 
PANELq 
QV 

RDLIN 

R1PRNT 

SHAPEG 

SHAPEL 

VERFN 

XI  NT 

VIBPOP 
VIBTMP 
UP PATH* 

/BNDRIES/ 
/CNTRL  / 
/PATHD  / 

FSCATM* 

/BNDRY  / 
/CONSTN/ 
/DEAMT  / 
/HMULS  / 
/PARMTR/ 
/ZOUTP  / 
(ATMCON) 
ATMPTH* 
PACK 
WATVAP 
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Obtains  Refractive  Path  and  Calculates  Absorber 


Atmospheric  Path  Driver 
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> 


FSCGEO 

ALAYER 

AMERCE 

ANDEX 

AUTLAY 

EXPINT 

FDBETA 

FINDSH 

FNDHNM 

HALFWD 

RADREF 

REDUCE 

RFPATH 

SCALHT 

MDLATM* 

/MDLATM/ 
/TRAC  / 
(MCATMB) 
CONVRT 
DEFALT 
NSMDL 
RDUNIT 


LOWTRN* 

/CNSTNS/ 
/LCRD1  / 
/LCRD2  / 
/LCRD3  / 
/LCRD4  / 
/MART  / 
/MDLZ  / 
/MODEL  / 
/RAIN  / 
/TITl  / 
/USRDTA/ 
j  /ZVSALY/ 

(TITLE  ) 

AERNSM 

/MDATA  / 
(MARDTA) 
(MDTA  ) 
(PRFDTA) 
AERPRF 
MARINE 
STDMDL 

CIRRUS 
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Sets  up  selected  Atmosphere 


LOWTRAN  Driver  for  Aerosols 


■» 


iV 

fv 

r 

y  4 
v, 

V 

V 

V, 

h  u 


1 
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EXAB1N*  Loads  the  Aerosol  Extinction,  Absorption  and 

Asymmetry  Coefficients  for  the  desired  Model 
/EXTD  /*  Stores  all  the  Aerosol  Attenuation  Coefficients 
(EXTDTA)*  Contains  all  the  Aerosol  Attenuation  Coefficients 
AB 

AITK 

DEBYE 

DOP 

GAMFOG 

GMRAIN 

INDX 


GEO 

/PARMLT/ 
/RFRPTH/ 
ANDEX 
EXPINT 
FDBETL 
FILL 
FINDSL 
FNDHML 
GEOINL 
LOLAYR 
PACKL 
RAD REF 
REDUCL 
RFPATL 
SCALHT 


NEWMDL 


TRANS*  Outputs  Aerosol  Optical  Properties  for  Path 

AB 

AEREXT*  Interpolates  Aerosol  Attenuation  Coeffients  for 

values  at  Wavenumber  v 

AITK 

DEBYE 

DOP 

GAMFOG 

GMRAIN 

INDX 

TNRAIN 


VSA 

PATH 

TLEMIS* 

EMOUT 


Merges  Multiple  Scattered  and  Single  Scattered 
Radiance  for  the  Total  Path 
Writes  Merged  Radiance  and  Transmittance  Results 


V. 


•'.V' 


.V' 
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XLAYMS* 


XMERGE* 


ABIN1T 

I 

ABStNT 


ABSMRG 


ADARSL 


EMINIT* 


AERF* 

BBFN 

EMIN* 

E3* 

FLUXLP* 

FLXADD* 

GETEXT* 

MSIN* 

MSOUT* 

RADFN 


GETEXT* 

RADFN 

XMERG1 


ABSOUT 


ABSOUT 


ABSOUT 


AERF* 

GETEXT* 


AERF* 

BBFN 

EMIN* 

EMOUT 

E3* 

FLXADD* 

MSIN* 

MSOUT* 


Contois  Layer  by  Layer  Calculation  for  Contributing 
Scattering  Layers  Above  Desired  Path 
Evaluates  Extinction,  Scattering,  and  Asymmetry 
Computes  Black  Body  Function 
Reads  Optical  Depths  and  Calculates  Tgas 
Computes  Exponential  Integ.al  E3 
Merges  Fluxes  from  Previous  and  Current  Layer 
Calculates  and  does  Adding  of  Fluxes  for  Current 
Layer 

Reads  Extinction,  Absorption,  and  Asymmetry 
Reads  Merged  Fluxes  from  Previous  Layer 
Writes  Merged  Fluxes  for  Current  Layer 
Computes  Radiation  Function 

Driver  for  Merging  Options 

Reads  Extinction,  Absorption,  and  Asymmetry 

Computes  Radiation  Function 

Driver  for  Merging  Routines  using  RAD1NT 

Initializes  Optical  Depth  Calculation 
Writes  Merged  Optical  Depth  Results 

Merges  Optical  Depths  using  General  4-point 

Interpolation 

Writes  Merged  Optical  Depth  Results 

Merges  Optical  Depths  using  Lagrange  4-point 

Interpolation 

Writes  Merged  Optical  Depth  Results 

Adds  Absorption  and  Scattering  to  Common  Blocks 
for  use  in  a  Transmittance  Run 
Evaluates  Extinction,  Scattering,  and  Asymmetry 
Reads  Extinction,  Absorption,  and  Asymmetry 

Initializes  Radiance  Calculation  with  First  Layer 
and  Drives  Layer  Flux  Calculation 
Evaluates  Extinction,  Scattering,  and  Asymmetry 
Computes  Black  Body  Function 
Reads  Optical  Depths  and  Calculates  t 
Writes  merged  Radiance  and  Transmittance  Results 
Computes  Exponential  Integral  E3 
Calculates  and  does  Adding  of  Fluxes  for  Current 
Layer 

Reads  Merged  Fluxes  from  Previous  Layer 
Writes  Merged  Fluxes  for  Current  Layer 
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FLUXLP* 


RADINT 


RADMRG* 


AERF* 

BBFN 

EMIN* 

E3* 

FLXADD* 

MS1N* 

MSOUT* 


AERF* 

BBFN 

EMBND 

EMIN* 

EMOUT 


AERF* 

BBFN 

EMIN* 

EMOUT 

E3* 

FLXADD* 

MSIN* 

MSOUT* 


Drives  Layer  Flux  Calculation  without  Radiance 
Evaluates  Extinction,  Scattering,  and  Asymmetry 
Computes  Black.  Body  Function 
Reads  Optical  Depths  and  Calculates  Tgas 
Computes  Exponential  Integral  E3 
Calculates  and  does  Adding  of  Fluxes  for  Current 
Layer 

Reads  Merged  Fluxes  from  Previous  Layer 
Writes  Merged  Fluxes  for  Current  Layer 

Merges  Radiance  and  Transmittance  using  General 
4-point  Interpolation 

Evaluates  Extinction,  Scattering,  and  Asymmetry 

Computes  Black  Body  Function 

Evaluates  Boundary  Contribution  to  Radiance 

Reads  Optical  Depths  and  Calculates  Tgas 
Writes  Merged  Radiance  and  Transmittance  Results 

Merges  Radiance  and  Transmittance  using  Lagrange 
4-point  Interpolation,  and  also  Drives  Layer 
Flux  Calculation 

Evaluates  bxtinction,  Scattering,  and  Asymmetry 
Computes  Black  Body  Function 
Reads  Optical  Depths  and  Calculates  Tgas 
Writes  merged  Radiance  and  Transmittance  Results 
Computes  Exponential  Integral  E3 
Calculates  and  does  Adding  of  Fluxes  for  Current 
Layer 

Reads  Merged  Fluxes  from  Previous  Layer 
Writes  Merged  Fluxes  for  Current  Layer 


Table  3-3 

Segment  Load  File  for  FASC0DE2  with  Multiple  Scattering 


FASCODE 

TREE 

FASCOD2-(LAYRS .LASER, SCANFN , PLOTT , 

,FLTRFN ,TEST) 

FASC0D2 

GLOBAL 

CONSTS , FILHDR .MAIN , IFIL .LAMCHN , LASIV , ADRIVE , 

.MSACCT, 

(MSCONS 

FASC0D2 

* 

INCLUDE 

FASCOD2 .BUFIN.BUFOUT 

LAYRS 

TREE 

XLAYER-(OPPATHS,OPDPTHS,XMERGES, FLUXES, TLEMIS.XLAYMS) 

XLAYER 

GLOBAL 

LINHDR,XBLF,ABSORA,SCATTA,ASYMMA,RMRG,KLTSCT 

XLAYER 

* 

INCLUDE 

XLAYER, MSCOPY 

OPPATHS 

TREE 

OPP ATH-( FSCATMA , LOWT , PATH) 

OPPATH 

CLOBAL 

CMTRL .BNDRIES , PATHD-SAVE 

OPPATH 

* 

INCLUDE 

OPPATH 

FSCATMA 

TREE 

FSCATM-( FSCGEO .MDLATM) 

FSCATM 

INCLUDE 

FSCATM, ATMPTH.ATMCON , PACK.WATVAP 

FSCATM 

GLOBAL 

CONSTN , HMOLS .PARMTR , DEAMT .ZOUTP , BNDRY-SAVE 

FSCGEO 

INCLUDE 

RFP ATH , ALAYER , AUTLAY .EXPINT, FINDSH , SCALHT , ANDEX , 

.RADREF 

.HALFWD, FSCGEO , FDBETA, FNDHMN , REDUCE, AMERGE 

MDLATM 

INCLUDE 

MDLATM , NSMDL , MLATMB , CONVRT , RDUN IT , DEFALT 

MDLATM 

* 

GLOBAL 

MLATM.TRAC-SAVE 

FLUXES 

TREE 

FLXDWN-( SRCFCN ) 

FLXDWN 

* 

INCLUDE 

FLXDWN ,RADFN .GEXEXT , FLXADD ,MSIN ,MSOUT , E3 , AERF , BBFN 

LOWT 

TREE 

LOWTRN-(GEO , EXABIN .TRANS .AERNS , VSA .CIRRUS .NEWMDL) 

LOWTRN 

INCLUDE 

TITLE 

LOWTRN 

GLOBAL 

LCRD 1 .LCRD2 .LCRD3 , LCRD4 .MODEL .CNSTNS , 

.MDLZ.ZVSALY ,MART 

.USRDTA.TITL, RAIN-SAVE 

GEO 

INCLUDE 

GEO .GEOINL , REDUCL , EXPINT , RADREF , SCALHT , ANDEX , 

.RFPATL 

,  FDBETL ,  FNDHML ,  FINDSL ,  FILL  ,LOl,AYR,PACKL 

GEO 

GLOBAL 

P ARMLT , RFRPTH-S  AVE 

EXABIN 

INCLUDE 

EXABIN , EXTDTA , INDX . DEBYE , DOP , AB , CAMFOG , 

,  AITK.GMRALN 

EXABIN 

GLOBAL 

EXTD 

TRANS 

INCLUDE 

TRANS , AEREXT .TNRAIN , DEBYE , DOP , INDX ,AB , 

.GAMFOG 

,  AITK .GMRAIN 

AERNS 

TREE 

AERNSM 

AERNSM 

INCLUDE 

AERNSM , AERPRF , PRFDTA , MARINE .MARDTA .STDMDL ,MDTA 

AERNSM 

i k 

GLOBAL 

MDATA-SAVE 

PATH 

★ 

INCLUDE 

PATH 

XMERGES 

TREE 

XMERGE-(ABINIT .ABSMRG , ABSINT.EMINIT .RADMRG .RADINT , 

.ADARSL 

,  FLL’XLP ) 

XMERGE 

INCLUDE 

RADFN, XMERGE, GETEXT , XMERGI 

ADARSL 

INCLUDE 

ADARSL, GETEXT, AERF 
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* 

OPDPTHS  TREE  0PDPTH-(HIRAC1 ,LBLF4 ,CONTNW,NLTE) 

OPDPTH  GLOBAL  LBLF, ABSORB .SCATTR 

* 

H1RAC1  INCLUDE  HIRAC1 ,SHAPEL,SHAPEG , VO ICON ,RDLIN .CNVFNV .PANEL, MOLEC , 

,QV , ABSOUT , VERFN , Rl PRN  T , XINT , BMOLEC , BHIRAC , BSHAPEL 
HIRACI  GLOBAL  FNSH-SAVE 

* 

LBLF 4  INCLUDE  LBLF4 .RDLIN4 .CONVF4 ,MOLEC , QV , SHRINK, VOICON .BMOLEC 

* 

CONTNW  TREE  CONTNM-(SL296,SL260 ,FRN296 ,FRNC02 , 

,N2C0NT,O3CHAP,O3HHT0,O3HHTI .03HHT2 .02C0NT) 

CONTNM  INCLUDE  CONTNM.XINT 

03 CHAP  INCLUDE  03CHAP ,03CH 

03HHT0  INCLUDE  O3HHTO.O3HH0 

03HHT1  INCLUDE  03HHTI ,03HH1 

03HHT2  INCLUDE  03HHT2.03HH2 

02C0NT  INCLUDE  02C0NT.02C 

SL296  INCLUDE  SL296.BS296 

SL260  INCLUDE  SL260.BS260 

FRN296  INCLUDE  FRN296.BFH20 

FRNC02  INCLUDE  FRNC02.BFC02 

N2C0NT  INCLUDE  N2CONT.BN2 

* 

NLTE  TREE  NONLTE-(HIRACQ,VIBTMP .VIBPOP) 

NONLTE  GLOBAL  VBNLTE 

HIRACQ  INCLUDE  HIRACQ , SHAPEL .SHAPEG , VOICON ,RDLIN .CNVFNQ , PANELQ , 

, MOLEC , QV , ABSOUT , VERFN ,R1 PRNT ,XINT , BMOLEC , BHIRAQ , BSHAPEL 
HIRACQ  GLOBAL  FNSQ-SAVE 

* 

ABINIT  INCLUDE  ABINIT .ABSOUT 

* 

ABSMRG  INCLUDE  ABSMRG .ABSOUT 

ABSINT  INCLUDE  ABSINT .ABSOUT 

* 

EMINIT  INCLUDE  EMINIT .EMIN .EMOUT .MSIN .MSOUT.BBFN ,AERF .FLXADD ,E3 

* 

RADMRG  INCLUDE  RADMRG , EMIN .EMOUT , MSIN , MSOUT , BBFN ,AERF , FLXADD ,E3 

it 

RADINT  INCLUDE  RADINT, EMIN .EMOUT , BBFN ,AERF .EMBND 

* 

SRCFCN  INCLUDE  SRCFCN .BETABS .MSEMIS 

it 

TLEMIS  INCLUDE  TLEMIS, EMOUT 

* 

XLAYMS  INCLUDE  XLAYMS .GETEXT .RADFN , FLUXLP , MSIN .MSOUT , EMIN , EMOUT , 

.FLXADD, BBFN, AERF.E3 
* 

FLUXLP  INCLUDE  FLUXLP, FLXADD, MSIN, MSOUT, BBFN, AERF, EMIN, E3 

it 

SCANFN  INCLUDE  SCANFN , SHAPEG , RDSCAN , SHRKSC , SHAPET .CONVSC , PANLSC , 

, CNVREC , SINCSQ , CKPRNT , INTER 
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* 


FLTRFN 

it 

INCLUDE 

FLTRFN , RDSC AN , CNVFLT 

TEST 

if 

INCLUDE 

TEST.BTEST 

PLOTT 

PLTFAS 

PLTFAS 

.ENDPLT 

* 

TREE 

GLOBAL 

INCLUDE 

PLTFAS-(BBSCLE .HEADER .AXEST .FPLINE) 
PLTHDR.AXISXY.YCOM, POINTS, TITLOC, NAME 

PLTFAS .LINT ,EXPT .XNTLOG .TEMPFN .TENLOG.MNMX.FSCLIN 

AXEST 

TREE 

END 

AXES-(AXISL,AXLOG,AX2) 

FASCOD2 

Table  3-4 

File  Usage 


UNIT  NAMES: 

TAPES  -  INPUT 

TAPE6  -  OUTPUT 

TAPE10  -  Output  file  containing  gas  optical  depths  for  current  laver. 
Labeled  KFILE. 

TAPE11  -  Output  file.  For  FASCODE  runs  without  MR,  it  contains  merged 
results  from  previous  layer.  For  multiple  scattering  runs,  it 
contains  the  merged  total  radiance  and  transmittance  for  the 
path.  Label  alternates  between  LFILE  and  MFILE.  Always  labeled 
LFILE  for  final  output. 

TAPE12  -  Output  file.  For  FASCODE  runs  without  MS,  it  contains  the  final 
merged  results.  For  multiple  scattering  runs,  this  file 
contains  the  total  thermal  emission  radiance.  Label  alternates 
between  MFILE  and  LFILE.  Always  labeled  MFILE  for  final  output. 

TAPE19  -  Intermediate  output  file,  containing  aerosol  coefficients  for 
path  above  desired  path.  Only  used  if  upper  path  is  defined. 
Labeled  IMSFIL. 

TAPE20  -  Intermediate  output  file,  containing  aerosol  coefficients  for 
desired  path.  Labeled  IEXFIL. 

TAPE21  -  Intermediate  output  file,  containing  total  layer  flux  values. 
Label  alternates  between  IRFILE  and  LRFILE. 

TAPE22  -  Intermediate  output  file,  containing  composite  layer  flux 

values.  Label  alternates  between  JSFILE,  KSFILE,  and  MSFILE. 

TAPE23  -  Intermediate  output  file,  containing  composite  layer  flux 

values.  Label  alternates  between  JSFILE,  KSFILE,  and  MSFILF.. 

TAPE24  -  Intermediate  output  file,  containing  total  layer  flux  values. 
Label  alternates  between  ISFILE  and  LSFILE. 

TAPEZ1!  -  Intermediate  output  file,  containing  composite  layer  flux 

values.  Label  alternates  between  JSFILE,  KSFILE,  and  MSFILE. 

TAPE26  -  Intermediate  output  file.  This  file  holds  optical  properties 

for  upper  path  for  the  upward  adding  loop,  and  for  the  downward 
loop  it  contains  the  final  multiple  scattered  contribution  to 
the  radiance  for  the  FASCODE  path.  Initially  labeled  IEF1LE, 
labeled  KEFILF,  at  end  of  run. 

TAPE27  -  Intermediate  output  file,  containing  the  multiple  scattered 

contribution  to  the  radiance-  Label  alternates  between  IEFILF,, 
JEFILE,  and  KEFILE. 

TAPE28  -  Intermediate  output  file,  containing  the  multiple  scattered 

contribution  to  the  radiance.  Label  alternates  between  IEFILE, 
JEFILE,  and  KEFILE. 
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UNIT  NAMES 


KFILE 

LFILE 


MFILE 


1MSF1L 

1EXFIL 

ISFILE 

JSFILE 

KSFILE 


LSFILE 

MSFILE 


IEFILE 


JEF1LE 

KEF1LE 


Contains  gas  optical  depths  for  current  layer. 

Contains  results  of  previous  layer  for  FASCODE  runs  without  mul¬ 
tiple  scattering.  Contains  final  merged  results  for  multiple 
scattering  run. 

Contains  finai  merged  results  for  FASCODE  runs  without  multiple 
scattering.  Contains  total  single  scattered  results  for  multi¬ 
ple  scattering  runs. 

Contains  aerosol  coefficients  for  path  above  desired  path. 
Contains  aerosol  coefficients  for  desired  path. 

Contains  total  composite  fluxes  for  previous  layer. 

Contains  downward  composite  fluxes  for  current  layer. 

Contains  upward  composite  fluxes  for  previous  layers  during  up¬ 
ward  adding  loop.  Contains  downward  composite  fluxes  for  previ¬ 
ous  layer  during  downward  adding  loop. 

Contains  total  composite  fluxes  for  current  layer. 

Contains  upward  composite  fluxes  for  current  layer  during  upward 
adding  loop.  Contains  upward  composite  fluxes  for  all  layers 
during  downward  adding  loop. 

Contains  _  for  upper  path  during  upward  aduing  loop.  Con¬ 
tains  multiple  scattering  contribution  to  the  radiance  from  the 
current  layer  during  downward  adding  loop. 

Contains  multiple  scattering  contribution  to  the  radiance  from 
the  previous  layers. 

Contains  total  multiple  scattering  contribution  to  the  path  ra¬ 
diance. 


3.1.5  Notes  on  the  Operation  of  FASC0DE2  with  Multiple  Scattering 

In  order  to  use  FASC0DE2  with  multiple  scattering,  it  would  be  helpful 
to  understand  a  little  about  how  and  why  it  does  the  calculation  the  way  it 
does.  Supplemental  user  instructions  are  listed  in  Table  3-5.  The  first 
thing  that  the  user  will  note,  is  that  a  new  flag  appears  in  Card  1.2.  This 
flag  Is  for  turning  off  and  on  the  multiple  scattering  sections  of  the  pro¬ 
gram.  As  outlined  in  the  user  instructions,  various  other  parameters  must 
also  be  set  so  that  a)  aerosols  are  included,  b)  radiance  is  calculated,  c) 
FSCATM  is  called,  and  d)  normal  merging  is  done. 

The  next  card,  is  where  an  understanding  of  how  the  calculation  is  per¬ 
formed  can  be  very  useful.  Card  1.2a  is  a  new  card  which  is  only  read  If 
IMS=1.  This  card  includes  four  quantities,  a)  TGRND,  which  is  the  temperature 
of  the  ground  (at  0.0  km),  b)  SEMIS,  which  is  the  emissivity  of  the  surface 
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(at  0.0  km),  c)  KMINMS,  the  minimum  height  for  multiple  scattering  (default  H 
0.0  km),  and  d)  HMAXMS,  the  maximum  height  for  multiple  scattering  (default  • 
15.0  km).  HMINMS  and  HMAXMS  determine  the  bounds  that  are  placed  on  the  mul¬ 
tiple  scattering  calculation.  HMINMS  determines  the  height  at  which  multiple 
scattering  will  start.  HMAXMS  determines  the  height  above  which  multiple 
scattering  paths  are  not  calculated.  If  the  maximum  of  the  desired  path  is 
below  HMAXMS,  the  program  will  define  a  path  from  the  maximum  of  the  desired 
path  to  HMAXMS.  If  the  maximum  of  the  desired  path  is  above  HMAXMS,  then 
HMAXMS  is  ignored,  but  the  entire  path  is  calculated  using  multiple  scatter¬ 
ing.  This  means  that  HMAXMS  only  checks  to  see  if  the  contribution  from  above 
should  be  added,  but  does  not  cut  off  the  multiple  scattering  calculation  at 
that  height.  For  the  case  where  the  minimum  of  the  desired  path  is  above 
HMAXMS,  then  no  multiple  scattering  is  performed  for  that  path. 

For  example,  if  the  user  defined  a  slant  path  from  4.0  to  10.0  km,  with 
HMINMS  =  0.0  km  and  HMAXMS  =  15.0  km,  the  program  would  first  define  a  new 
path  from  10. 0  to  15.0  km,  to  obtain  che  downward  flux  contribution.  Second¬ 
ly,  a  path  from  0.0  to  4.0  km  would  be  defined  to  obtain  the  contribution  from 
the  layers  below  4  km.  Then  lastly  the  program  would  define  the  path  from  4.0 
to  10.0  km  and  calculate  the  desired  radiances,  adding  in  any  contributions 
obtained  from  the  first  two  calculations.  It  should  therefore  be  clear  that  a 
simple  path  from  4.0  to  10.0  km  using  multiple  scattering,  may  actually  in¬ 
volve  the  execution  of  FASC0DE  for  THREE  separate  paths.  This  therefore  can 
increase  the  run  time  proportionally. 

This  example  just  described,  is  the  worst  case  scenario.  Another  case 
could  consider  a  path  from  10  to  30  km,  with  HMINMS  and  HMAXMS  as  before.  The 
program  would  first  define  a  path  from  0.0  to  10.0  km  to  obtain  the  contri¬ 
bution  from  below,  and  would  then  do  the  path  from  10  to  30  kra  using  multiple 
scattering.  Yet  another  possibility,  is  to  define  a  path  from  0.0  to  10.0  km, 
with  HMINMS  and  HMAXMS  as  before.  The  program  would  first  do  a  path  from  10 
to  15  km  to  obtain  the  contribution  from  above,  and  then  do  the  path  from  0.C 
to  10.0  km.  The  final  possible  case,  is  to  define  a  path  from  0.0  to  20.0  km 
with  HMINMS  and  HMAXMS  still  set  to  0.0  and  15.0  km  respectively.  This  would 
result  in  only  one  path  from  0.0  to  20.0  km  calculating  the  multiply  scattered 
radiances . 
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One  thing  the  user  should  note,  is  that  in  the  preceeding  examples  there 
was  no  mention  of  angles,  or  tangent  paths.  That  is  because  this  methodology 
is  true  without  regard  to  the  angle  for  slant  paths,  and  is  true  for  all  tan¬ 
gent  paths.  For  example,  if  you  want  to  consider  a  tangent  path  starting  at 
25  km  with  a  minimum  at  10  km,  and  a  range  of  1000  km,  leaving  HM1NMS  and 
HMAXMS  as  previously  defined.  The  program  would  first  do  the  slant  path  from 
0,0  to  10  km  to  obtain  the  contribution  from  below,  and  would  then  do  the  de¬ 
sired  tangent  path. 

The  following  items  should  be  noted  by  the  user: 

1)  When  doing  a  multiple  scattering  run,  the  input  file  is  rewound  so 
it  is  not  possible  to  string  multiple  scattering  jobs  together. 

2)  Additional  paths  which  are  defined  by  the  multiple  scattering 
parameters  will  increase  run  time  proportionaly . 

3)  Some  coding  notes: 

a)  The  coding  which  sets  and  determines  the  multiple 
scattering  paths  is  contained  in  ATMPTH. 

b)  The  two  routines,  XLAYMS  &  FLUXLP  were  written  for  the 
multiple  scattering  paths  which  don't  require  the 
calculation  of  radiances. 

A  test  case  for  FASC0DE  along  with  the  input  and  output  files  is 
described  in  Appendix  B. 
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Table  3-5 

Supplemental  User  Instructions  for  Using  FASC0DE2 
with  Multiple  Scattering 


The  following  user  instruction  cards  are  not  a  complete  FASC0DE2  user 
deck,  but  only  present  the  changes  to  existing  cards  and/or  the  addition  of 
new  cards  required  for  Multiple  Scattering  runs. 


1)  CARD  1.2  now  appears  as  follows: 


1HIRAC, 

ILBLF4, 

1CNTNM, 

1AERSL, 

1EM1T, 

ISCAN, 

IFILTR, 

I PLOT, 

5. 

10, 

15, 

20, 

25, 

30, 

35, 

40, 

4X.I1, 

4X.I1, 

4X.I1, 

4X.I1, 

4X.11, 

4X.I1, 

4X.11, 

4X.11, 

ITEST, 

IATM, 

IMRG, 

I  LAS, 

IMS, 

1RAD, 

MPTS , 

NPTS 

45, 

50, 

55, 

60, 

65, 

70, 

71-75, 

76-80 

4X.I1, 

4X.11, 

4X.11, 

4X.11, 

4X.I1, 

4X.I1, 

15, 

15 

where  each  of  the  parameters  is  as  defined  in  the  FASC0DE2  user  instructions, 
excluding  the  new  flag,  IMS  which  is  used  for  turning  on  the  multiple 
scattering  sections  of  FASC0DE2, 


IMS  (0,1)  flag  for  multiple  scattering 

■  0  normal  FASCODE  run,  no  scattering 
-  1  FASCODE  run  with  multiple  scattering 

NOTE:  the  following  must  also  be  true  for  multiple  scattering: 

a)  1AERSL  =  1  or  7 

b)  1EMIT  =  1 

c)  IATM  =  1 

d)  IMRG  =*  0 

2)  CARD  1.2a  (this  card  directly  follows  Card  1.2) 


ste  sfc  *  if 

TGRND , 

1-10, 

E10.3. 


Required  if  IMS  =  1,  otherwise  omit. 


SEMIS, 

KMINMS , 

HMAXMS 

11-20, 

21-30, 

31-40 

E10.3. 

E10.3. 

El  0 . 3 

Table  3-5  (Cont.) 


TGRND  Temperature  of  surface  (degrees  Kelvin) 
(no  default  value) 

SEMIS  Emissivity  of  surface 
(no  default  value) 

HMINMS  Minimum  height  for  Multiple  Scattering 
(default  =  0.0  Km) 

HMAXMS  Maximum  height  for  Multiple  Scattering 
(default  =  15.0  km) 


3)  CARD  3.6  (optional  card  for  LOWTRN  subroutine)  now  appears  as  follows: 

(1  to  16)  (DUMMY, EXTC(1, I), ABSC( l ,1) ,ASYM( 1 , I), (1*1,47)) 

Format  (3(F6.2,2F7 .5.F6.4) ) 

where  EXTC,  and  ABSC  are  as  defined  in  FASC0DE2  user  instructions,  and  ASYM 
is  the  Aerosol  Asymmetry  Factor  which  is  required  for  multiple  scattering 
calculations . 

NOTE:  When  multiple  scattering  runs  are  performed,  subsequent  runs  cannot  be 
done  using  only  one  Input  deck.  This  is  due  to  the  rewinding  of  the 
input  deck  during  multiple  scattering  execution.  It  is  therefore 
required  that  each  multiple  scattering  run,  necessitates  its  own  input 
deck. 
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3 . 2  LOWTRAN  Implementation 

3.2.1  Overview 

In  this  section,  we  describe  the  technical  details  of  incorporating 
into  LOWTRAN  6  the  multiple  scattering  treatment  together  with  the  k-distri- 
bution  method.  The  procedure  and  the  new  code  structure  are  also  presented. 
The  existing  LOWTRAN  6  program  structure  is  shown  in  Figure  3-3  (Figure  27  of 
Kneizys  et  al.,  1983).  For  simplicity  of  implementation,  we  have  added  only 
two  new  subroutines,  MSRAD  and  FLXADD,  to  handle  the  k-distribution  and 
multiple  scattering.  The  subroutine  FLXADD  is  called  from  MSRAD.  The 
functions  of  MSRAD  and  FLXADD  are  shown  in  Figure  3-4,  while  the  individual 
elements  are  described  below. 

3.2.2  k-Dlstrlbutlon  Method 

Because  of  the  complicated  molecular  absorption  band  structure  of  the 
gases,  a  rigorous  frequency  integration  would  require  a  line-by-line 
integration  which  clearly  is  very  time  consuming  and  thus  not  acceptable  for 
LOWTRAN.  Our  approach  Lo  this  problem  is  the  application  of  the  absorption 
coefficient  (k)  distribution  method  (cf,  Wang  and  Ryan,  1983;  Stephens, 

1984).  In  a  homogeneous  gas  1^_ er,  the  k-distribution  function  is  formally 
related  to  the  mean  transmission  function  T^Cu), 

T  (u)  =  -r—  /  e-ku  d v  =  /  f(k)e-kudk 

Av  A v  J0 

(3.34) 

=  / o  e~kudg  *  l  e  kiU  Agi 
1  =  1 

where  Av  is  the  narrow  repeated  interval  (20  cm-1  in  LOWTRAN)  and  u  is  the  gas 
amount.  The  f(k)  for  a  given  gas  at  a  specified  Av  is  the  probability  density 
function  such  that  f(k)dk  is  the  fraction  of  the  frequency  interval  for  which 
the  absorption  coefficient  is  between  k  and  k  +  dk.  Eqn.  (3.34)  reveals  that 
the  transmission  depends  on  the  distribution  of  k-values  within  Av  but  not  on 
the  ordering  of  the  values.  The  cumulative  k-distribution  function  is  g(k), 
while  (k^,  Ag^)  are  the  discreet  sets  of  values  to  approximate  the  integral. 
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Figure  3-3.  LOVTTRAN  6  Main  Program  Structure 
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TRANS 

•  Absorption  parameters  for  gases 


Aerosol  scattering  and  extinction  parameters 


MS  RAD 

•  Aerosol  extinction  and  scattering  optical  thickness 

•  Rayleigh  scattering  optical  thickness 

•  Asymmetry  factors 

•  Radiance  adding  using  source  functions  from  FLXADD 


FLXADD 

•  Molecular  optical  thickness:  H20,  uniformly  mixed  gases,  ozone 
©  Multiple  scattering  two-streara  flux  parameterization  for  solar 

and  thermal  for  each  homogeneous  layer. 

•  Adding  method  for  fluxes  at  each  layer  for  inhomogeneous  atmosphere 
©  Solar  and  thermal  source  functions 


Figure  3-4.  Functions  of  subroutine  TRANS,  MSRAD,  and  FLXADD  used  in  the 
LOWTRAN  implementation  of  multiple  scattering. 
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For  th.'  . umraation  in  Eqi  .  (3.34)  to  be  useful  with  LOWTRAN,  values  for  k 
and  Ag  must  be  found  to  fit  the  LOWTRAN  transmission  data  for  water 
vapor/uniformly  mixed  gases  (H20+)  and  ozone  (O3).  The  transmission  for  water 
vapor/unif ormly  mixed  gases  and  ozone  may  be  expressed  as 


Ti 

t2 


M  . 

T(H20+)  =  l  e  iUl  Agi 
i=l 
N 

=  T(03)  =  l  e'kJU2  Ag j . 
3=1 


(3.35) 

(3.36) 


The  problem  of  fitting  LOWTRAN  transmission  data  as  an  exponential  sum  has 
been  handled  successfully  by  Wiscombe  and  Evans  (cf.  Wiscombe  and  Evans, 

1977),  using  10  k  and  Ag  values  for  H20+  and  six  for  O3.  Their  k  values  and 
associated  Ag  values  are  listed  in  Table  3-6.  Plots  of  the  percent  transmis¬ 
sion  error  of  their  fit  relative  to  the  LOWTRAN  data  for  both  H20+  and  O3  as  a 
function  of  loglO  of  the  effective  optical  depth  are  shown  in  Figures  3-5  and 
3-6. 


Making  use  of  the  Wiscombe  and  Evans  exponential  sum  fit  of  the  LOWTRAN 
transmission  data  yields  values  of  k^,  kj  ,  and  Ag^,  Agj  for  M  =  10,  N  ■  6. 
However,  if  M  ^  N  and  Ag^  ^  Agj ,  it  is  necessary  to  perform  multiple 
scattering  calculations  separately  for  water  vapor/uniformly  mixed  gases 
( H2 0+)  and  ozone  (O3 )  using  two  probability  (k)  loops.  This  problem  is 
avoided  by  choosing  either  the  10  H20+  or  6  O3  Ag^  values  and  then  expressing 
the  transmission  of  both  H20+  and  O3  as  the  sum  of  either  10  or  6  terms.  Run 
time  efficiency  concerns  suggest  that  both  transmissions  be  expressed  as  the 
sum  of  six  exponentials,  using  the  Ag.  values  for  the  ozone  fit  shown  in  Table 
3-6.  Thus,  in  the  LOWTRAN  implementation,  the  transmission  function  for  water 
vapor  is  refit  to  the  Ag  values  for  ozone  for  each  path.  The  transmission  as 
a  function  of  probability  is  shown  in  Figure  3-7.  (Figure  is  illustrative 
only  —  actual  values  are  not  plotted.) 
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Fig.  3-5  Percent  error  of  transmittance  fit  for  water  vapor  and  uniformly 

mixed  gases  relative  to  LOWTRAN  values.  Transmittance  fit  is  based 
on  Wiscombe  and  Evans,  1977 
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Figure  3-7.  Transmission  vs.  cumulative  probability  for  Wiscombe  and  Evans 

(1977)  fit  to  LOWTRAN  H20+  and  O3  transmission  values.  (Figure  does 
not  depict  actual  transmission  values.) 


As  an  example,  the  refit  H2O  transmission  for  the  first  of  six  components  may 
be  expressed  as 


T1  1  21  [Tll  ASll  +  T12  ASl2  +  T13  A2l3  +  T14^821  “  813)]/  2 1  • 

.  6 

where  Tu  *  e  iUl  ,  i  =  1,  10,  so  that  T(H20+)  =  XI  =  l  Agj_.  The  total 

molecular  transmission  for  H2Q+  and  O3  for  a  particular  layer  for  some  small 

frequency  interval  is  now 


T1  •  t2 


[  £  igi] 

i-1 


[  £  e-kl'“2  Jgl]. 

i-1 


(3.38) 


Assuming  no  overlap  of  gas  absorption,  the  many  cross  terms  may  be  eliminated, 
resulting  in 


T1 


•  T 


2 


l  e‘(kiUl  +  kiU2)  Agl. 

i-1 


(3.39) 


Figure  3-8  compares  the  product  of  the  water  vapor  and  ozone  obtained 
from:  (1)  the  intrinsic  LOWTRAN  band  model,  (2)  the  k-distribution  fit 

including  overlap  terms  (Eqn.  (3.38)),  and  (3)  the  k-distribution  fit  assuming 
no  overlap  terms  (Eqn.  (3.39)),  i.e.,  eliminating  cross  terms.  This 
calculation  is  for  the  near  surface  layer  of  a  US  standard  atmosphere  run  and 
covers  the  spectral  range  from  about  6.7  to  20  pm.  It  can  be  seen  that  there 
is  some  difference  between  the  calculated  transmission  but  very  little  of 
practical  consequence  for  the  multiple  scattering  calculation. 

The  individual  components  of  the  combined  gas  (H20+  +  O3)  molecular 
absorption  transmission  can  thus  be  expressed  by 


T i  “  T j £  •  T 2  £ ,  i  -  1 ,  2 ,  . . .  6 .  (3. AO) 

The  jth  component  of  the  optical  thickness  of  a  given  layer  due  tu  molecular 
absorption  is 


i  =  -  In 
u 


<V- 


(3.41) 


Tj  is  added  to  the  molecular  continuum,  aerosol  extinction,  and  molecular 
scattering  optical  thickness  to  determine  the  Lotal  optical  thickness  lor  each 
k  value  and  layer  in  FLXADD. 


TRANS 


It  is  important  to  note  that  the  k-distribution  transmittances  are  em¬ 
ployed  only  in  the  determination  of  the  multiply  scattered  source  function. 
The  adding  of  radiances  based  on  the  source  function  depends  on  the  original 
LOWTRAN  transmissions.  Thus,  errors  in  the  combined  k-distribution/streara 
approximation  approach  for  multiple  scattering  are  not  propagated  beyond  the 
approximate  multiple  scattering  approach. 


Table  3-6 

The  Set  of  k.£  and  Values  Used  to  Approximate  the 
LOWTRAN  6  Transmission  Data 


i 

ki 

Agi 

Water 

1 

Vapor/Uniformly  Mixed  Gases 

1.718754  x  10" 2 

1.233690  x 

10"4 

2 

2.490766  x  10"! 

4.614758  x 

10"  3 

3 

2.328222 

2.404693  x 

10"  2 

4 

7.744676  x  10"* 

2.772890  x 

10“  2 

5 

2.229876  x  10"* 

5.409872  x 

10"  2 

6 

4.438688  x  10"2 

2.119191  x 

\ 

o 

.H 

7 

1.081388  x  10"2 

2.776904  x 

10"  1 

8 

3.054533  x  10"3 

3.089278  x 

10"  1 

9 

1.258653  x  10“3 

9.072606  x 

10"  2 

10 

0.000000 

1.240893  x 

10"4 

Ozone 

1 

4.736517  x  10*1 

7.395116  x 

10"2 

2 

3.823806  x  10“2 

5.557539  x 

10"  1 

3 

7.424869  x  10"3 

2.476504  x 

10"  1 

4 

1.785435  x  10~3 

8.688930  x 

10"2 

5 

4.336295  x  10"4 

3.565867  x 

IQ"2 

6 

0.000000 

1.385517  x 

10"4 

II, 


'Xv 

u 
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3.2,3  Inhomogeneous  Atmosphere 

For  inhomogeneous  atmospheres,  we  adopt  the  same  scaling  approximation 
used  in  LOWTRAN  (see  Eqn.  2.21): 


k.(P,G)  =  k.(P  ,9  )  — /  0  /0 

1  1  O  0  r  O 

o 

We  use  the  same  adding  method  discussed  in  Section  3.1.2  to  calculate 
the  thermal  radiation  flux  in  an  inhomogeneous  atmosphere.  Basically,  the 
parameters  used  are  F+,  F~,  T,  and  R,  as  presented  in  Section  3.1.2. 


3.2.4  Stream  Approximation,  Source  Function,  and  Radiance  Calculation 

Again,  we  use  the  same  procedures  discussed  in  Section  3.1.3.  to 
calculate  the  radiance  from  the  source  function  and  stream  approximation. 

The  multiple  scattered  contribution  to  the  source  function  is 
approximated  by  (Eqn.  2.14,  2.2,  2.3,  2.10,  2.11). 

Jms(t  ,±u  ,$  )  -  -jp- —  j  lx(r  )  1 1-6  (y  ) )  +  F+(t)  B  (  u  ) }  •  (2.14) 

The  total  source  function  is: 


•J  (t  ,u  ,<}i  )  "  J0(r  ,u  A  )  +  JMS^T  •$) 


(2.2) 


where : 


•J0(t  ,u  »4> )  = 


+  [  1-!0o(t  )]  B [G (t  )] 


(2.3) 


The  single  scattered  contribution  to  .1  above  is  taken  directly  from  the 
existing  LOWTRAN  single  scattering  coding.  (Actually  the  summed  multiply 
scattered  radiance  is  added  to  the  summed  singly  scattered  radiance.) 

One  essential  difference  between  the  multiple  scattering  radiance 
implementation  and  the  previous  single  scattering  version  is  the  treatment  ot 
surface  reflection.  Through  the  surface  boundary  condition,  surface 
reflection  affects  the  flux  profile  and  hence  source  function  throughout  the 
atmosphere.  From  the  radiance  solutions: 


Bv? 


'■*  -  V“  .  V  “ »  V  . ’«  **»  ■  *»’  * .  *■_  \ 

.  J*  _  .  J  _  J*.  A..*  _ 


I(t  ,ii  ,iji )  -  {—  [iry  Fe”1  Uo  +  /  /  I(T*,-y,<j>)ydyd$] 


+  (1-r)  B[0(t*)}  exp-(T*-x)/y  +  f  J(t,y,$)  e"(t“T)/lJ  ^ 
-  /  J(t,y,*)  e"^-t)/y  dt 


(2.10 


(2.11 


It  can  be  seen  that  surface  reflection  can  increase  both  upward  and  downward 
radiances  through  the  source  function.  For  upward  radiances,  there  is  a  more 
drastic  difference  between  the  standard  LOWTRAN  treatment  and  the  implemented 
multiple  scattering  version.  This  concerns  reflection  of  downward  scattered 
radiance  from  the  surface  and  back  to  a  downward  looking  observer  (the  second 
term  in  the  brackets  in  Equation  (2.10)).  This  contribution  is  not  included 
in  standard  LOWTRAN  although  the  reflected  attenuated  direct  solar  term  is 
(i.e.,  the  first  terra  in  Equation  (2.10)). 

3.2.5  Description  of  L0WTRAN6  with  Multiple  Scatterin 


The  addition  of  multiple  scattering  to  L0WTRAN6  has  been  made  with  a 
minimum  of  changes  to  the  source  code  and  input  deck.  Only  one  new  input 
parameter,  1MULT,  has  been  added.  IMULT,  read  as  the  laat  input  on  card  1,  is 
set  to  zero  for  the  original  LOWTRAN6  with  single  scattering  and  one  for  the 
AER  L0WTRAN6  with  multiple  scattering.  If  aerosols  are  not  included  (clear 
sky  case),  the  original  L0WTRAN6  is  executed,  independent  of  the  value  of 
IMULT. 

Source  code  changes  were  made  primarily  in  the  main  program  LWTRN6  and  in 
subroutines  TRANS,  with  some  additions  to  EXABIN,  AEREXT,  EXTDTA,  and  RFPATH . 
Subroutines  MSRAD,  FLXADD,  and  ALEVEL  and  functions  BETABS  and  PLANCK  were 
added . 

1 ne  main  program  has  been  altered  to  allow  the  path  geometry  calculations 
(.ailed  l rum  subroutine  GEO)  to  be  accessed  twice,  once  for  the  original 
LOWTKAN  si  gle  scatter  path  and  once  for  the  entire  atmosphere,  required  for 
imputation  oi  the  multiply  scattered  radiance.  The  reason  for  the  second 
path  geometry  call  is  that  the  thermal  and  solar  source  functions  used  to 
determine  the  multiply  scattered  radiance  contribution  are  functions  of  the 
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upward  and  downward  flux  at  each  layer.  The  multiply  scattered  flux  at  a 
given  layer  will,  in  part,  be  determined  by  radiation  from  all  atmospheric 
layers,  even  those  above  or  below  the  layers  between  the  observer  and  target. 
To  calculate  the  fluxes  at  each  layer,  it  is  therefore  necessary  to  add  the 
flux  contributions  from  the  surface  up  to  space,  and  back  down  to  the  surface 
again.  This  adding  procedure,  executed  in  the  new  subroutine  FLXADD,  is  dis¬ 
cussed  in  detail  in  Sections  3. 1.2.1  and  3. 1.2.2. 

A  flag  (variable  ITEST)  has  been  incorporated  into  RFPATH  for  the  purpose 
of  isolating  the  refracted  viewing  path  zenith  angle  for  each  layer,  as  op¬ 
posed  to  the  solar  zenith  angle,  calculated  for  the  solar  cases  in  a  later 
pass  through  RFPATH. 

A  number  of  changes  were  incorporated  into  subroutine  TRANS.  After  cal¬ 
culating  the  cumulative  path  parameters  used  in  the  original  L0WTRAN6,  it 
calls  MSRAD,  where  the  aerosol  scattering  and  extinction,  H2O  continuum  and 
Rayleigh  scattering  optical  thickness  for  each  layer  are  calculated.  The 
asymmetry  factor  for  each  layer  Is  also  calculated  based  on  model  extinction 
and  absorption  data  added  to  EXAB1N,  AEREXT,  and  EXTDTA,  and  aerosol  effective 
absorber  amounts  for  the  four  vertically  spaced  aerosol  regions.  MSRAD  calls 
subroutine  FLXADD. 

Subroutine  FLXADD  calculates  the  six  degraded  k  components  of  the  H2O/ 
uniformly  mixed  gas  transmission  for  each  layer,  as  well  as  the  molecular  ab¬ 
sorption  optical  thickness.  It  is  then  added  to  the  continuum,  aerosol  ex¬ 
tinction,  and  molecular  scattering  optical  thickness  from  MSRAD  to  provide  the 
total  optical  thickness  for  each  layer  and  k-value,  as  well  as  the  correspond¬ 
ing  single  scatcer  albedo.  The  diffuse  flux  contribution  for  each  isolated 
layer  and  k-value  is  computed  from  the  applicable  two-stream  approximation, 
either  solar  or  thermal,  and  combined  in  the  flux  adding  routine  to  determine 
the  total  upward  and  downward  flux  for  each  layer  for  each  k  value.  The  dif¬ 
fuse  solar  downward  flux  is  summed  over  k  value  and  returned  to  evaluate  the 
surface  reflected  downward  diffuse  solar  flux.  Function  PLANCK  returns  the 
black  body  radiance  for  a  particular  wavenumber  and  temperature  in  units  of 
Wctn  2  sterad  ^cm  1  ,  while  subroutine  ALEVEL  returns  the  layer  number  corre¬ 
sponding  to  a  particular  height,  the  top  layer  of  the  atmosphere  being  layer 
1,  Knowing  the  upward  and  downward  fluxes  as  well  as  the  backscatter  param¬ 
eter  8  returned  from  function  BETABS,  the  multiple  scattered  source  function 
can  be  determined. 
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Returning  to  MSRAD,  the  radiance  is  summed  over  the  viewing  path  (between 
HI  and  H2)  and  k-values  to  provide  the  multiply  scattered  diffuse  radiance 
contribution.  For  a  downward  looking  (upward  radiance)  calculation,  the 
multiply  scattered  radiance  at  the  upper  boundary  may  be  expressed  by 


6  n,  —  1 

I+(nt)  =  l  {(1-T^)  S(k, nt )  +  b£  (T{+1  -  TL)  S  (k.i  +  l)j  Agk  +  Is  (3.39) 
k=l  i=nt 

where  is  the  bottom  layer,  nt  the  top  layer,  S(k,i)  the  source  function  for 

a  particular  k  value  and  layer  i.  is  the  transmission  from  the  upper 

boundary  (H2)  through  layer  i  along  the  viewing  path.  The  total  upward 

radiance  also  includes  I  ,  a  direct  solar  reflection  surface  contribution  or 

s 

surface  thermal  emission  term,  which  for  the  new  AER  LOWTRAN  includes  the 
contribution  due  to  the  reflection  of  the  single  and  multiple  scattered  solar 
radiance.  For  a  thermal  case, 

I  =  B  T  e  (3. AO) 

S  8  t 

where  Bg  is  the  black  body  function  for  a  given  surface  temperature  and 
frequency,  Tt  is  the  total  atmospheric  transmissivity,  and  e  is  the  surface 
emissivity.  For  the  old  LOWTRAN  (single  scatter)  solar  case, 

I  *  a  y  T  S/tt  ,  (3.41) 

s  sot 

where  ag  5s  the  surface  albedo,  u0  the  cosine  of  solar  zenith  angle  at  the 
surface,  and  S  the  direct  solar  intensity  at  the  surface.  With  multiple 
scattering  (IMULT=1),  the  single  and  multiple  scatter  reflection  terms  are 
added,  so  that 


I  «  a  y  TS/ir  +  a  FT  /*,  (3.42) 

s  sot  s  s  t 

where  FT  is  the  downward  diffuse  flux  at  the  surface. 

S 

For  an  upward  looking  (downward  radiance)  case,  the  multiply  scattered 
radiance  is  expressed  by 


nt+1 


l~(n  )  =  l  {(1  -  T  )  S(k,n  )  +  £  (T  -  T  )  S  (k,  1-1)}  Ag.  ,  (3.43) 

o  11  n,  b  ,  II  -1  K. 

k=»l  b  i^n, 

b 


where  T^  is  the  transmission  from  the  lower  boundary  (HI)  through  layer  i 
along  the  viewing  path.  At  this  point,  the  L0WTRAN6  path  parameters  are 
reloaded,  and  TRANS  is  executed  once  more  for  the  single  scatter  L0WTRAN6 
case.  The  single  scatter  diffuse  radiance  is  added  to  the  multiple  scatter 
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term  and  the  boundary  radiance  contributions  (if  any)  to  yield  the  total 
thermal  radiance,  while  the  solar/lunar  radiance  is  the  sum  of  the  single  and 
multiple  scatter  terms  plus  a  surface  reflection  terra  where  applicable.  MSRAD 
and  FLXADD  are  not  called  for  the  L0WTRAN6  single  scatter  case. 


Figures  3-9  and  3-10  are  flow  charts  of  subroutines  TRANS,  MSRAD,  and 


FLXADD . 


LOWTRAN  6 
Set  up  atmo¬ 
spheric  model 


Subroutine  TRANS 
Being  frequency  loop 
v  “  v  +  Av 

1"*"  A.  X 

-  T . . 

Molecular  absorption 
coefficients 


Layer  Loop 

Load  appropriate  absorber 
amounts  for  thermal/solar  use 


Scattering?]-^- 
i  I  MULT  -  l?l>-- 


|  Calculate  H2Q/unif orraly  mixed  gases, 
O3 ,  HjO  continuum  and  aerosol 
cumulative  transittance 


Solar  ? 


lerraal  radiance  addins 


Subroutine  GEO:  _ 

Calculate  and  store  absorber  amts 
for  single  scatter  optical  path 

Subroutine  GEO: 

Calculate  and  store  absorber 
amounts  for  complete  atmosphere 
for  multiple  scattering  use _ 

Subroutine  MSRAD: 

H2O  continuum  and  aerosol  optical 
thickness  for  each  layer 


Asymmetry  Factor  Calculation  j _ 

_ x  .Subroutine 

_  '  1  FLXADD  __ 

— ~~1  Start  k  loop 

— 1  Start  layer  loop 

Multiply  Scattered  Radiance 
adding  for  thermal/solar 

I  End  of  k  and  layer  loops  I 


Single  scatter  , 
_V  solar  radiance  | 
^  adding  I 


II _  _  _  yes_ 
ISet  IMULT  =  0 


[MULT  = 


1  -  -  -  -|End  Frequency  Loot 


Figure  3-9.  Flowchart  of  Subroutines  TRANS  and  MSRAD. 


From  MSRAD 


^  j  Subroutine  FLXADD  1 


Begin  K  loop  | 


Calculate  ground  layer  flux 


Begin  layer  loop  for  upward  adding  | 


Calculate  tmol(N ,k) ,Ttot(N,k)  =  Tmol(N,k)  +  Taer(N) 
Single  scatter  albedo  W(N,k)  =  Tg^gCt  (N)  /Ttot  (N,k) 


Two-stream  approximations:  thermal  and  solar 
Layer  and  composite  transmissivity,  reflectivity, 
Flux  calculations:  Ft,  Ft  _  _ 


_  _  |  Begin  downward  adding  layer  loo 


I  1  Complete  calculation  of  Ft  and  Ft  for  each  layer 


End  downward  loo 


Thermal  and  solar  source  functions  | 


.-I  End  of  k  loop] 


.  —|  Return  to  MS  RAD 


Figure  3-10.  Flowchart  of  Subroutine  FLXADD. 
3.2.6  Notes  on  the  Operation  of  L0WTRAN6  with  Multiple  Scatterin 


As  previously  discussed,  the  only  change  to  the  LOVfTRAN  input  card 
images  occurs  in  card  1  with  the  addition  of  the  input  parameter  IMULT  used  to 
turn  multiple  scattering  on  or  off.  Card  1  is  now  read  as  follows: 

MODEL , 1TY PE , IEMSCT ,M1 , M2 , M3 , 1M , NOPRNT ,TB0UND , SALB , IMULT 
15  ,  15  ,  15  ,15,15,15,15,  15  ,F10.3  ,F10.3,  15 

As  previously  mentioned,  setting  IMULT  =  1  initiates  multiple  scattering, 
while  IMULT  =  0  provides  single  scattering  or  clear  sky  results. 
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Unlike  the  original  L0WTRAN6,  the  new  LOWTRAN&  requires  a  complete 
atmospheric  profile  (surface  to  space)  to  properly  calculate  the  multiple 
scattering  contribution  to  the  radiance.  For  a  calculation  involving  a  path 
from  0-1  km,  for  example,  it  is  still  necessary  to  choose  the  desired 
stratospheric  as  well  as  higher  level  aerosol  distribution  set  by  1VULCN  on 
card  2.  For  a  3-20  km  path,  the  boundary  layer  aerosols,  set  using  IHAZE,  on 
card  2  or  user-defined  via  card  2D,  will  now  affect  the  calculated  radiance. 
The  original  LOWTRAN  performing  single  scattering  required  aerosol  information 
only  along  the  path  between  observer  and  target. 

As  currently  set  up,  L0WTRAN6  will  not  do  multiple  Rayleigh  ccattering 
for  IHAZE  =  0  but  will  instead  default  to  the  original  LOWTRAN6  single 
scattering  calculation.  To  perform  a  multiple  Rayleigh  scattering  calculation 
without  boundary  layer  aerosols,  choose  IHAZE  >  0  and  set  the  visibility  (VIS 
on  card  3)  to  a  large  value  such  as  300  km  to  effectively  eliminate  the  effect 
of  these  aerosols. 

When  L0WTRAN6  is  run  in  its  original  f*.  (single  scattering  only, 

IMULT  =  0),  the  output  for  a  solar  case  includes  the  reflected  attenuated 
direct  solar  radiance  under  the  heading  GROUND  REFLECTED.  However,  in  many 
cases,  the  ground  reflected  diffuse  radiance  due  to  single  and  multiple 
scattering  is  o.  the  same  order  of  magnitude  as  the  direct  reflected 
radiance.  The  output  for  new  L0WTRAN6  (IMULT  *  1)  has  therefore  been  altered 
to  include  both  the  total  (diffuse  single  and  multiple  scattered  and  direct) 
and  direct  reflection  terms.  Since  the  original  single  scatter  LOWTRAN  does 
not  calculate  downward  radiance  at  the  surface  for  an  upward  radiance  case, 
the  single  scattered  diffuse  reflection  term  cannot  currently  be  calculated. 

Running  LOWTRANo  with  multiple  scattering  results  in  a  several-fold 
increase  in  CPU  time.  The  atmospheric  profile  and  path  being  used  may  not 
produce  a  significant  multiple  scattering  contribution,  possibly  resulting  in 
wasted  computer  time.  A  new  feature  has  therefore  been  added  to  the  solar 
output.  Roth  the  total  (single  and  multiple)  and  single  scattered  radiance 
are  now  included  under  the  heading  "PATH.  SCATTERED  RADIANCE."  If  in  doubt, 
the  user  may  now  run  a  multiple  scatter  LOWTRAN  run  to  determine  the  relative 
aize(s)  of  the  single  and  multiple  scatter  contributions.  If  multiple 
scattering  is  negligible,  future  calculations  may  be  made  with  the  more 
economical  single  scatter  (IMULT  =  0)  option. 

Test  cases  for  LOWTRAN  are  provided  in  Appendix  C. 


4.  RESULTS  OF  INTERCOMP ARISON  OF  IMPLEMENTED  MULTIPLE  SCATTERING 
APPROACH  WITH  EXACT  CALCULATIONS 


4.1  Comparison  of  FASCODE  Multiple  Scattering  Results  to  Exact  Calculations 

The  detailed  approach  described  in  the  previous  sections  was  employed  to 
develop  a  multiple  scattering  version  of  FASCODE  based  on  the  stream 
approximation  and  adding  method.  In  order  to  verify  the  operation  of  the 
code,  a  number  of  test  cases  were  run.  To  validate  the  results,  these  test 
cases  were  also  run  with  the  discrete  ordinate  code  (DOM)  for  multiple 
scattering.  Additionally,  the  existing  nonscattering  version  of  FASCODE  was 
run  for  comparison.  The  optic.  properties  for  theae  test  cases  and 
corresponding  radiance  results  are  given  in  Table  4-l(a)-(c).  Three  test 
cases  were  selected.  Case  1  is  based  on  a  FASCODE  run  at  900  cm-1  with  a 
rural  haze  and  5  km  visibility.  This  resulted  In  an  optical  depth  of  0.357 
and  a  single  scattering  albedo  of  0.126.  With  this  small  albedo,  the  effects 
of  scattering  are  expected  to  be  small  and  results  for  all  three  methods 
should  be  comparable.  In  order  to  enhance  the  scattering  effects  and 
emphasize  the  differences  between  approaches,  the  number  density  of  aerosol 
was  increased  to  nonphysical  amounts  to  yield  single  scattering  albedos  in  the 
vicinity  of  0.5  (case  II  )  and  0.9  (case  III).  Note  that  the  FASCODE  models 
include  the  effect  of  atmospheric  refraction  while  the  discrete  ordinate 
method  code  is  strictly  a  plane  parallel  model.  Comparisons  of  these  results 
for  zenith  angles  in  the  domain  from  65  to  115  degrees  should  be  made  with 
this  in  mind. 

Results  are  plotted  in  Figures  4-1  to  4-6.  Upward  and  downward  emergent 
radiances  for  Case  I  are  shown  in  Figures  4-1  and  4-2,  respectively.  Here 
FASCODE,  FASCDDEMS,  and  the  DOM  are  shown  to  give  approximately  the  same  mag¬ 
nitude  and  angular  distribution  of  radiance.  A  slight  tendency  for  FASCODE  to 
underestimate  the  upward  radiance  and  overestimate  the  downward  radiance  can 
be  noted.  At  a  higher  optical  depth  and  single  scattering  albedo  (case  II), 
this  effect  becomes  obvious  (Figures  4-3  and  4-4,  respectively.)  Here  it  can 
be  seen  that  FASCODEMS  does  a  clearly  superior  job  at  representing  the  angular 
radiance  distribution.  For  a  situation  where  multiple  scattering  dominates 
over  gas  absorption  (case  III,  Figures  4-5  and  4-6),  the  deficiency  of  the 
absorption  only  FASCODE  model  is  obvious. 
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Table  4-1 


Scattering  Comparisons  for  FAS CODE  Multiple  Scattering  Runs 

(WAVELENGTH  -  000  CM-1  ,  RADIANCES  -  ERG  SEC-1  CM-2  (CM-l)-l  STR-1) 
(a)  CASE  I.  SMAL'j  SCATTERING 


TAU=0.357,  OMEGA0=0 .126,  ASYMMETRY  FACTOR-0.6985 


ZENITH  ANGLE 
(DEGREES) 

AFGL  FASCODE 
(NO  MS) 

AER  FASCODE 
(MS) 

DISCRETE 
ORDINATE  METHOD 
(MS) 

UPWARD 

RADIANCES 

0.000 

101.29 

102.55 

106.24  t 

8.349 

10C 

102.42 

103.07 

19.165 

1 '  . 

101 .84 

102.45 

30.045 

V 

100.69 

101.23 

40.939 

.29 

98.75 

99.18 

51.839 

94.02 

95.58 

95.81 

62.741 

88.62 

90.23 

90.09 

73.644  * 

79.17 

80.66 

79.61 

84.548  * 

68.23 

68.29 

60.96 

DOWNWARD 

RADIANCES 

95.452  * 

65.16 

63.44 

68.51 

106.356  * 

46.79 

44.37 

46.72 

117.259 

35.15 

32.85 

34.09 

128.161 

28.48 

26.36 

27.12 

139.061 

24.43 

22.48 

23.01 

149.955 

21.92 

20.11 

20.49 

160.835 

20.41 

18.70 

19.00 

171.651 

19.64 

17.98 

18.24 

180.000 

19.47 

17.81 

14.69  t 

tThis  point  extrapolated  in  discrete  ordinate  method, 

*FASC0DE  path  is  refracted  while  discrete  ordinate  method  is  plane  parallel. 
These  zenith  angles  are  affected. 
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(2)  CASE  II.  INTERMEDIATE  SCATTERING 

TAU-0.555,  OMEGA0-0.498,  ASYMMETRY  FACTOR-0.6985 


ZENITH  ANGLE 
(DEGREES) 

AFGL  FASCODE 
(NO  MS) 

AER  FASCODE 
(MS) 

DISCRETE 
ORDINATE  METHOE 
(MS) 

UPWARD 

RADIANCES 

0.000 

94.73 

101.62 

103.52  t 

8.349 

94,55 

101.48 

101.43 

19.165 

93.77 

100.85 

100.69 

30.045 

92.25 

99.63 

99.21 

40.939 

89.79 

97.53 

96.68 

51.839 

86.00 

94.02 

92.50 

62.741 

80.24 

88.21 

85.45 

73.644  * 

71.93 

78.66 

73.33 

84.548  * 

68.01 

68.06 

55.37 

DOWNWARD  RADIANCES 

95.452  * 

66.48 

59.53 

70.14 

106.356  * 

56.04 

44.76 

50.95 

117.259 

45.63 

34  .02 

37.25 

128.161 

38.49 

27.40 

29.14 

139.061 

33.78 

23.35 

24.28 

149.955 

30.72 

20.89 

21.32 

160.835 

28.84 

19.46 

19.58 

171.651 

27.87 

18.71 

18.70 

180.000 

27.65 

18.55 

16.80  t 
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(c)  CASE  III.  LARGE  SCATTERING 


TAU ”0.308, 

OMEGAO^O .897 ,  ASYMMETRY 

FACTOR-O.6985 

ZENITH  ANGLE 
(DEGREES) 

AFGL  FASCODE 
(NO  MS) 

AER  FASCODE 
(MS) 

DISCRETE 
ORDINATE  METHOD 
(MS) 

UPWARD 

RADIANCES 

0.000 

103.31 

112.24 

119.78  t 

8.349 

102.97 

112.19 

113.19 

19.165 

102.46 

111.96 

112.81 

30.045 

101.46 

111.49 

112.03 

40.939 

99.78 

110.53 

110.58 

51.839 

97.01 

108.57 

107.87 

62.741 

92.28 

104.51 

102.49 

73,644  * 

83.37 

95.13 

90.55 

84.548  * 

68.85 

68.79 

62.28 

DOWNWARD  RADIANCES 

95.452  * 

63.04 

46.54 

57.61 

106.356  * 

41.52 

21.81 

28.65 

117.259 

30.60 

12.71 

16.29 

128.161 

24.75 

8.72 

10.68 

139.061 

21.33 

6.78 

7.84 

149.955 

19.24 

5.83 

6.31 

160.835 

18.00 

5.37 

5.48 

171,651 

17.37 

5.14 

5.09 

180.000 

17.23 

5.09 

- t 

In  order  to  gain  some  insight  into  the  differences  between  the  FASCODE 
and  FASCODEMS  results,  a  radiance  budget  has  been  computed  for  one  upward  and 
one  downward  zenith  angle  corresponding  to  case  II.  These  radiance  budgets 
are  presented  in  Table  4-2  and  provide  the  magnitude  of  contributions  to  the 
total  radiance  from  surface  emission,  atmospheric  emission,  and  the  scattering 
source  function,  as  applicable.  Radiance  values  are  rounded  for  simplicity. 

Table  4-2 

Case  II  Radiance  Budgets  [unattenuated  surface  |missipn  corresponds  to 
120  radiance  units;  units  =  erg  s  cm  str  (cm  )  ]. 


FASCODE 

FASCODEMS 

DOM 

(a)  Upward  Radiance 

(zenith  angle  8,3 

degrees ) 

Surface  emission 

65 

65 

Atmospheric  emission 

30 

16 

Scattering  source 

20 

Total 

95 

101 

101 

(b)  Downward  Radiance 

(zenith  angle  171 

.7  degrees) 

Atmospheric  emission 

28 

12 

Scattering  source 

- 

8 

Total 

28 

20 

19 

Examination  of  Table  4-2  indicates  that  FASCODE  over  estimates  atmospher¬ 
ic  emission  in  the  presence  of  scatterers  since  it  assumes  that  scattering 
extinction  is  equivalent  to  absorption.  In  this  case  radiation  absorbed  is 
re-emitted.  While  scatterers  may  absorb  and  re-emit,  some  fraction  of  the 
extinction  due  to  their  presence  (i.e.,  that  given  by  their  single  scattering 
albedo)  is  sitaply  redirected.  In  the  case  of  the  upward  intensity  in  Table 
4-2,  forward  scattering  of  some  of  the  surface  emission  lost  due  to  extinction 
more  than  compensates  for  the  decreased  atmospheric  emission.  Thus  the  non¬ 
scattering  version  of  FASCODE  underestimates  the  radiance.  For  the  downward 
radiance  in  Table  4-2,  the  small  scattering  source  does  not  compensate  for  the 
decrease  in  atmospheric  emision  due  to  the  treatment  of  scattering.  In  this 
case  FASCODE  overestimates  the  actual  radiance. 
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Figure  4-5. 


Comparison  of  FASCODE,  FASCODEMS,  and  DOM  for  upward  radiance 
[units:  ergs  sec  1  cm  2  str  1  (cm 
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4.2  Comparison  of  LOWTRAN  Multiple  Scattering  Results  Co  Exact:  Calculations 


4.2.1  Solar  Multiple  Scattering 

In  order  to  verify  the  correct  operation  of  the  implemented  LOWTRAN 
multiple  scattering  treatment  for  both  solar  and  thermal  regimes,  comparisons 
were  made  to  exact  results. 

In  the  case  of  solar  multiple  scattering,  exact  results  were  obtained  for 
a  variety  of  cases  during  the  trade-off  analysis  summarized  in  Appendix  A. 
Exact  solutions  to  the  radiative  transfer  equation  for  solar  multiple 
scattering  were  obtained  using  the  Gauss-Seidel  iterative  method  based  on  a 
code  by  Dave  (1972).  This  algorithm  evaluates  the  fluxes  and  radiances  for 
inhomogeneous  atmospheres  with  arbitrary  vertical  distributions  of 
anisotropically  scattering  aerosol  overlying  a  Lambert  reflecting  surface. 
Emergent  fluxes  and  radiances  from  this  code  were  compared  to  those  obtained 
from  LOWTRAN  for  three  values  of  total  optical  thickness  (0.25,  0.50,  1.00) 
and  two  values  of  surface  albedo  (0.0,  0.4). 

A  number  of  special  modifications  of  LOWTRAN  were  necessary  to 
accommodate  the  comparison.  These  included  reading  in  the  same  aerosol 
extinction  and  absorption  coefficients,  asymmetry  factor,  and  phase  function 
as  were  used  for  the  Dave  model  runs,  and  forcing  LOWTRAN  to  use  the 
quantities  with  the  same  vertical  distribution  as  the  Dave  code.  (Ordinarily, 
user  supplied  aerosol  data  read  into  LOWTRAN  will  only  be  used  in  the  lowest 
2  km  of  the  atmosphere.)  These  changes  made  the  comparison  as  close  as 
possible  to  a  direct  one.  The  LOWTRAN  code  was  run  at  0.55  ym  with  the  sun  at 
a  zenith  angle  of  60°. 

The  first  test  of  the  LOWTRAN  implementation  is  whether  the  adding  method 
is  calculating  the  correct  fluxes  for  use  in  the  stream  approximation  of  the 
multiply  scattered  source  function.  Table  4.3  compares  the  emergent  fluxes 
calculated  by  the  exact  code  and  LOWTRAN.  For  simplicity  all  fluxes  and 
radiances  plotted  here  have  been  normalized  to  correspond  to  it  units  of 
incident  solar  irradiance.  To  obtain  the  appropriate  engineering  units 
(watts/cm2  era  1,  and  watts/cm2  str  cm  1,  for  flux  and  radiance,  respectively), 
the  normalized  values  can  be  multiplied  by  1.719  x  10  s.  As  can  be  seen  from 
the  results,  the  solar  two  stream  approximation  and  flux  adding  procedure  in 
LOWTRAN  is  reproducing  the  exact  values  quite  well  throughout  the  range  of 


optical  thicknesses  and  surface  albedos  examined.  The.  errors,  on  the  order  of 
a  few  percent,  are  to  be  expected  from  the  two  stream  approximation. 

Figures  4-9  through  4-12  illustrate  the  comparison  between  the  exact 
radiance  calculation  as  a  function  of  path  zenith  angle  (theta)  and  those  ob¬ 
tained  from  LOWTRAN  using  either  the  solar  single  scattering  option  (LOWTRAN) 
or  the  new  multiple  scattering  option  (LOWTRAN  MS).  Plotted  are  the  normal¬ 
ized  emergent  upward  and  downward  radiances  for  an  optical  depth  of  0.5  and 
surface  albedos  of  0.0  and  0.4,  The  plots  extend  to  a  maximum  zenith  angle  of 
60°  since  at  larger  angles  the  effects  of  refraction  in  the  LOWTRAN  code  make 
it  difficult  to  compare  to  the  plane  paraLlel  exact  results.  In  all  cases  it 
can  be  seen  that  the  introduction  of  the  multiply  scattered  contribution  to 
path  radiance  in  LOWTRAN  MS  has  considerably  improved  the  simulation  as  com¬ 
pared  to  that  obtained  from  the  exact  code.  The  improvement  is  considerable, 
especially  for  upward  radiance  when  there  is  surface  reflection  (see  Figure  4- 
11).  This  is  because  although  the  single  scatter  version  of  LOWTRAN  included 
the  surface  reflection  of  the  attenuated  direct  solar  beam,  it  contained  no 
provision  to  treat  the  reflected  downward  scattered  radiance.  This  contribu¬ 
tion  (which  is  the  second  term  In  the  brackets  in  Equation  (2.10))  may  be  con¬ 
siderable  when  multiple  scattering  is  a  factor  and  is  included  in  Che  LOWTRAN 
MS  version  since  the  downward  scattered  flux  is  calculated  as  a  result  of  the 
adding  method. 

The  relative  accuracies  as  a  funtion  of  path  zenith  angle  expressed  as 
percent  errors  obtained  in  the  comparisan  between  LOWTRAN  MS  and  the  exact 
calculations  for  all  optical  depths  evaluated  are  summarized  in  Figures  4-13 
and  4-14,  for  surface  albedos  of  0.0  and  0.4,  respectively.  In  general,  the 
solar  multiple  scattering  approach  implemented  within  LOWTRAN  underestimates 
radiance  by  10  or  20  percent.  These  accuracies  are  consistent  with  those  ob¬ 
tained  off  line  in  the  trade-off  analysis  (see  Appendix  A,  Table  A-5.) 


4.2.2  Thermal  Multiple  Scattering 


Froper  operation  of  the  code  for  thermal  multiple,  scattering  was  veri¬ 
fied  by  comparing  LOWTRAN  MS  to  corresponding  FASCODE  MS  results.  FASCODE  MS 
has  been  compared  to  exact  multiple  scattering  using  the  discrete  ordinate 
method  in  Section  4.1.  The  case  chosen  was  a  vertical  path  from  the  surface 
to  20  km  using  a  U.S.  Standard  atmosphere  (MODEL  b)  and  a  rural  aerosol  model 


with  5  km  visual  range  (1HAZE  2).  The  spectral  domain  for  both  runs  was  850- 
1.150  cm  1 .  Figures  4-1 5a,  b  illustrate  this  comparison  for  upward  and  down¬ 
ward  radiance,  respectively.  The  FASCODE  MS  results  have  been  degraded  to 
5  cm  -  for  consistency  with  the  LOWTRAN  MS  values.  In  general  both  models 
agree  quite  well  within  this  spectral  region.  The  slight  discrepancy  between 
degraded  FASCODE  and  LOWTRAN  spectra  (especially  for  downward  radiance)  in 
this  region  is  not  due  to  multiple  scattering  effects,  but  rather  to  differ¬ 
ences  in  the  corresponding  clear  air  spectra.  This  is  illustrated  by  Figure 
4-16  where  all  aerosols  have  been  removed  from  the  calculation. 

One  problem  with  the  previous  treatment  of  thermal  scattering  within 
LOWTRAN  was  the  failure  to  provide  a  source  function  to  introduce  multiply 
scattered  radiance  contributions  along  the  observed  path.  As  a  consequence  of 
this  ad  hoc  approach,  LOWTRAN  seriously  underestimated  path  radiance  for  long 
paths  such  as  those  near  the  horizon  where  multiple  scattering  could  contri¬ 
bute  significantly  (Ben-Shalom  et  al.,  19  This  deficiency  is  corre.cted 

witn  the  implementation  of  the  multiple  scattering  approach  descibed  here.  To 
illustrate,  radiances  have  been  simulated  for  a  ground  based  observer  scanning 
the  sky  between  8.0  and  13.5  pm.  Figures  4-17,4-18,  and  4-19  illustrate  the 
calculated  results,  for  observing  zenith,  10°  above  the  horizon,  and  the  hori¬ 
zon,  using  the  standard  LOWTRAN,  the  Ben-Shalom  modification  (Ben-Shalom  et 
al,  1980),  and  LOWTRAN  MS.  For  comparison,  the  black  body  function  of  the 
lowest  atmospheric  layer  is  also  plotted.  The  calculation  has  assumed  a  mid¬ 
latitude  summer  atmosphere  and  5  km  visual  range. 

As  can  be  seen  in  Figure  4-17,  there  is  very  little  difference  in  the 
results  among  the  three  approaches  when  observing  zenith  when  multiple  scat¬ 
tering  effects  are  minimized.  For  larger  air  mass  factors,  however,  the  three 
approaches  diverge,  with  the  LOWTRAN  MS  result  bounded  below  by  the  standard 
LOWTRAN  and  above  by  the  Ben-Shalom  modified  version  (see  Figure  4-18).  This 
is  understandable  since  standard  LOWTRAN  removes  all  scattered  radiance  from 
the  path  radiance  and  the  Ben-Shalom  modification  assumes  that  all  scattered 
radiance  is  returned  to  the  path  radiance  (i.e,,  the  scattering  is  conserva¬ 
tive).  Actually  some  of  the  scattered  radiance  is  returned  to  the  path  radi¬ 
ance  through  the  scattering  source  function.  The  fraction  returned  defends  on 
the  single  scattering  albedo  and  througn  the  scattering  phase  function,  on  the 
viewing  geometry.  At  the  horizon  (Figure  4-19),  the  standard  LOWTRAN  result 
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significantly  underestimates  the  radiance  level  expected  qualitatively,  i.e., 
that  approximating  the  emission  of  the  lowest  atmospheric  layer. 

In  Figure  4-20a  are  plotted  LOWTRAN  MS  results  for  the  case  described 
above  with  the  observer  path  at  various  elevation  angles.  Tor  comparison, 
Figure  4-2o-  reproduces  a  set  of  radiance  measurements  taken  with  qualitative¬ 
ly  similar  conditions  (Bell  et  al,1960).  Notf  the  wavelength  range  of 

triese  measurements  is  much  broader  than  that  illustrated  in  Figure  4-20a  and 
the  vertical  scales  are  different.  Additionally,  the  ambient  temperature  when 
the  measurements  were  taken  is  about  5  K  higher  than  that  of  the  surface  layer 
of  the  raidlatitude  summer  model  used  in  the  LOWTRAN  MS  simulation.  In  spite 
of  these  differences,  the  simulated  and  measured  radiances  exhibit  similar 
magnitudes  and  dependence  on  path  elevation  angle. 
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Table  4-3 

Comparison  of  solar  multiply  scattered  emergent  fluxes 
(normalized  to  r  units  of  incident  irradiance)  from 
LOWTRAN  FLXAOD  subroutine  (L)  and  exact  calculation  (E). 


Upward  Flux 

F+  (n  -  0) 

-a  m  -m  m.m.  auiaia.  aia  m  •  k  »  mm~m  m  m.-  »*»  m.  m 

downward  Flux 

F  (it  *  n*) 

7T  * 

r  =  0.0  r  =  0.4 

L  /  E  I.  /  E  ! 

r  =  0 . 0  r  “  0 . 4 

L  /  E  _  L  /  E 

0 . 2 r) 

.188  /  .191  .611  /  .634 

.323  /  .316  .367  /  .371 

0 . 50 

.244  /  .256  .562  /  .584 

.503  /  .489  .536  /  .544 

1 , 0  u 

.316  /  .324  .492  /  .515 

1 

.599  /  .567  .605  /  .615 
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L0WTRAN6  downward  radiance  (solar,  .55  pm)  with  and  without 
multiple  scattering  compared  to  exact  (Dave)  results.  Sola 
zenith  angle  is  60°. 
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4-13,  Percent  error  of  LOWTRAN6  multiple  scatter  calculations 

relative  to  exact  (Dave)  resulta  for  optical  depths  between 
0.25  and  1.00  with  a  surface  reflectance  of  0,0. 
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^  Figure  4-14.  Percent  error  of  LOOT  RAN  6  multiple  scatter  calculations 

rv  relative  to  exact  (Dave)  res  Its  for  optical  depths  between 

0.25  and  1.00  with  a  surface  reflectance  of  0.4. 
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-1 5a.  Comparisons  of  LOWTRAN  and  FASCODE  upward  thermal  radiance 
(with  multiple  scattering)  between  850  and  1150  cm  . 
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4-1 5b.  Comparison  of  LOWTRaN  and  FASCODE  downward  thermal  radiance 
(with  multiple  scattering)  between  850  and  1150  cm  . 
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Figure  4-20a • 


Figure  4-20b. 


Downward  tiir-rrep.l  radiance  results  from 
LOW THAN  wirh  multiple  scattering  (8.0 
to  13.5  um).  Viewing  path  elevation  angle 
varies  between  0  to  00“ .  Also  included 
is  the  black,  bony  radiance  of  the  lowest 


atmospheric  layer  (  C). 


Downward  radiance  measurements  for 
various  viewing  angles  under  acmes 
pheric  conditions  qualitatively 
similar  to  tnose  of  Figure  <* .  20a 
(from  Bell  et  al.,  1960;. 


SUMMARY  AND  CONCLUSIONS 


D . 


5 . I  Summary 

The  goal  of  this  research  has  been  to  provide  a  suitable  parameterization 
of  multiple  scattering  for  implementation  within  the  LOWTRAN  and  FASCODE 
atmospheric  transmittance/radiance  raoaels  developed  at  AFGL.  Due  to  the  com¬ 
prehensive  spectral  domains  treated  by  these  models  it  was  desired  to  develop 
a  uniform  aporoach  applicable  to  spectral  regions  from  the  ultraviolet  to  the 
millimeter  wave.  Furthermore  for  consistency,  analogous  approaches  were 
desired  for  implementation  in  LOWTRAN  and  FASCODE.  The  specific  nature  of  the 
chosen  approach  was  also  constrained  by  the  requirements  of  desired  accuracy 
and  computational  efficiency  of  these  models  and  additionally  by  their 
inherent  code  structure.  Code  structure  requirements  were  particularly 
limiting  with  respect  to  the  FASCODE  model. 

Aftec  reviewing  the  relevant  radiative  transfer  literature  two  potential 
candidate  multiple  scattering  parameter izations  were  selected  for  testing. 
These  were  a  stream  approximation  and  a  successive  order  of  scattering  ap¬ 
proach.  Extensive  incercoraparison  tests  were  conducted  to  compare  the  accur¬ 
acy  and  efficiency  of  these  candidate  approaches  to  exact  numerical  treatments 
in  the  solar,  thermal,  and  mictowave/mi 1 limeter  wave  scattering  regimes. 

Based  on  the  results  of  these  i ntercomparisons  it  was  decided  that  the  stream 
approximation  would  provide  the  required  degree  of  accuracy  and  the  least  com¬ 
putational  level  of  effort. 

Two  additional  considerations  were  factors  in  formulating  appropriate 
approaches  to  introduce  multiple  scattering  Into  the  LOWTRAN  and  FASCODE 
models.  First,  in  the  case  of  LOWTRAN  it  was  necessary  to  provide  an 
interface  between  the  molecular  Lsorber  band  model  treatment  and  tne  multiple 
scattering  parameterization.  To  accomplish  this  task  the  ^-distribution 
method  was  introduced.  The  k-d l at r l but i on  method  effectively  decouples  the 
m:lii;>le  scattering  f  i  om  the  LOWTRAN  treatment  of  i  ngray  gas  absorption  so 
that  t  tie  spectral  integration  can  oe  done  properly. 

Secondly  ,  alt  lough  multiple  scaltejing  approaches  generally  begin  with  a 
knowledge  ot  the  monuchroraat i c  optical  ptopefties  of  the  entire  path  of  inter¬ 
est,  both  K.V,  a:vi  Y  A.-NlOUil  pertorm  l  i  ,nsmi  t  lance  ami  radiance  ca  1  cu  la  t  i  on  s 

on  a  1  -L>,--li  ..i  ‘m  ,  i  s .  Thus  ,  in  t  a.-sv  models,  the  properties  of  the 


overall  path  are  not  known  until  the  end  of  the  calculation.  To  further 
complicate  matters,  the  spectral  resolution  of  each  FASCODE  layer  is 
determined  by  the  appropriate  pressure  dependent  average  line  width.  Thus, 
for  a  given  desired  spectral  domain,  the  resolution  of  each  FASCODE  layer 
increases  with  decreasing  pressure  and  monochromatic  optical  data  on  the 
entire  desired  scattering  path  is  never  available.  These  characteristics  of 
the  FASCODE  code  structure  led  to  the  adoption  of  a  layer  adding  approach  in 
order  to  accommodate  both  the  stream  approximation  and  the  layer-by-layer 
nature  of  the  calculation. 

Path  dependent  radiances  at  900  cm  1  were  calculated  using  the  FASCODE 
implementation  of  the  multiple  scattering  parameterization  (FASCODEMS),  the 
original  FASCODE,  and  the  discrete  ordinate  method.  Comparison  of  these 
results  indicates  that  .he  multiple  scattering  version  of  FASCODE  provides 
more  accurate  radiance  distribution  in  the  presence  of  aerosol  scattering. 

Our  parametric  studies  indicate  that  accuracies  of  IS  to  20  percent  should  be 
expected.  However,  this  increase  in  accuracy  incurs  an  additional  cost  in 
code  execution  time.  In  our  timing  test,  CPU  execution  times  for  a  multiple 
scattering  run  in  a  three-hundred  wave  number  region  (800-1100  cm”1),  took 
about  3-1/2  times  longer  than  a  standard  FASCODE  run.  Most  of  this  additional 
time  was  due  to  the  required  spectral  merging  of  stored  quantities  rather  than 
calculations  of  the  stream  approximation.  Additionally,  multiple  scattering 
runs  of  FASCODE  require  three  times  the  temporary  storage  of  regular  FASCODE 
runs.  The  program  size  is  increased  by  about  15  percent. 

To  insure  consistency,  LOWTRAN  and  FASCODE  multiple  scattering  results 
were  also  compared.  A  much  closer  agreement  is  found  between  AER  LOWTRAN  and 
FASCODE  (both  use  multiple  scattering)  than  those  between  L0WTRAN6  (uses 
single  scattering)  and  FASCODE. 

5.2  Recommendations 

This  report  has  described  recommendations  and  implementations  of  an 
approach  to  provide  the  LOWTRAN  and  FASCODE  models  with  a  unified  treatment  of 
multiple  scattering  consistent  with  described  levels  of  accuracy  and  their 
respective  code  structures.  Introduction  of  this  capability  has  a  number  of 
implications  which  require  additional  consideration.  These  form  the  basis  of 
our  recommendations. 


For  FASCODE,  the  question  of  solar  multiple  scattering  remains  an  open 
issue.  Certainly,  the  line  parameter  compilation  provides  the  capability  to 
exercise  FASCODE  well  below  2500  cm”1  where  solar  effects  can  influence  path 
radiance.  We  did  not  Introduce  solar  scattering  into  FASCODE  since  the 
necessary  solar  irradiance  and  single  scattering  routines  from  LOWTRAN  are  not 
yet  a  part  of  FASC0DE2 .  If  solar  multiple  scattering  is  desired  in  FASCODE, 
these  routines  and  the  capability  to  calculate  the  direct  solar  contribution 
and  the  scattering  phase  function  (i.e.,  the  solar  path  transmission)  should 
be  added  to  FASCODE,  then  a  suitable  coupling  to  the  stream  approximation/flux 
adding  code  segments  could  be  implemented. 

Likewise,  although  this  report  has  illustrated  the  applicability  of  the 
scattering  parameterization  to  the  microwave  and  millimeterwave  region  for 
scattering  by  precipitation,  the  FASCODE  model  does  not  provide  the  necessary 
optical  properties  for  such  calculations.  Attenuation  (i.e.,  extinction) 
coefficients  for  rain  are  available  in  FASCODE;  however,  information  on  single 
scattering  albedo  (or  alternatively,  the  scattering  coefficients)  and 
asymmetry  factor  are  not.  It  would  seem  a  simple  matter  to  provide  this 
data.  A  simple  parameterization  of  these  values  in  the  domain  between  lv.35 
and  231  GHz  (0.65  -  7.7  cm"1)  is  available  in  the  interim  (Isaacs  et  ai., 
1985). 

With  regard  to  LOWTRAN ,  we  recommend  that  the  k-values  should  be  part  of 
the  laboratory  experiments,  i.e.,  they  be  provided  directly  from  fitting  the 
laboratory  measurements  of  transmission  values.  The  present  setup  of  LOWTRAN 
can  easily  accommodate  the  updated  k-values  of  different  gases.  The  scaling 
approximation  to  treat  inhomogeneous  atmospheric  effects  needs  to  also  be 
revised,  perhaps  utilizing  the  correlated-k  techniques  published  by  Wang  and 


Finally,  it  is  worth  noting  that  although  the  choice  of  the  stream 
approximation  for  multiple  scattering  for  implementation  within  the  LOWTRAN 
and  FASCODE  models  was  not  an  arbitrary  one,  other  approaches  are  possible. 

One  criteria  which  the  user  may  wish  to  consider  is  the  trade-off  between 
accuracy  and  computational  level  of  effort.  As  implemented,  the  multiple 
scattering  enhancements  described  here  will  provide  emergent  radiances  with 
percent  errors  generally  in  the  range  of  10  to  30  percent  for  solar  scattering 
and  less  than  10  percent  for  thermal  scattering.  Internal  radiances,  i.e., 
those  evalauted  for  paths  within  the  atmosphere,  will  be  somewhat  less 
accurate.  These  accuracies  should  be  adequate  for  many  simulation 
requirements.  For  users  requiring  higher  accuracies,  the  code  changes  to 
LOWTRAN  and  FASCODE  to  accommodate  multiple  scattering  have  made  it  possible 
for  experienced  users  of  these  models  to  access  the  necessary  optical 
properties  profiles  for  input  to  numerical  multiple  scattering  algorithms. 

This  can  be  accomplished  by  further  modifying  the  LOWTRAN  and  FASCODE  models 
to  produce  output  files  of  the  profiles  of  optical  thickness,  single 
scattering  albedo,  asymmetry  factor,  and/or  angular  scattering  function  which 
are  used  in  the  set  up  of  the  current  stream  approximation  calculation. 

Other  possible  options  which  may  interest  knowledgeable  users  include 
output  of  emergent  fluxes  or  flux  profiles  (made  possible  by  the  inclusion  of 
the  two  stream/adding  step  done  in  the  stream  approximation)  and  substituting 
other  multiple  scattering  approximations  (such  as  higher  order  stream 
approximations)  within  the  code  to  obtain  more  accuracy.  One  application 
which  we  have  not  explored  is  the  applicability  of  the  implemented  code 
changes  to  evaluation  of  multiply  scattered  visible  and  infrared  radiances  by 
clouds.  A  variety  of  exciting  avenues  are  opened  with  the  inclusion  of  the 
multiple  scattering  capability  within  LOWTRAN  and  FASCODE. 
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APPENDIX  A 

REVIEW  AND  TRADE-OFF  ANALYSIS  OF  MULTIPLE  SCATTERING  PARAMETERIZATIONS 


A . I  Review  of  Multiple  Scattering  Parameterlzatlons 

Radiative  transfer  theory  provides  the  mathematical  description  of  the 
interaction  between  both  external  (e.g.,  incident  solar)  and  internal  (e.g., 
emitted  thermal)  sources  of  electromagnetic  radiation  and  the  optically  active 
atmospheric  constituents  of  gases  and  particulates  (e.g.,  aerosols,  clouds, 
and  precipitation).  The  solution  to  the  radiative  transfer  equation  for  a 
scattering  atmosphere  (Chandrasekhar,  1960;  Goody,  1964;  Liou,  1980)  generally 
requires  a  numerical  solution.  An  extensive  variety  of  exact  numerical 
methods  have  been  developed  to  solve  these  equations  over  the  years,  and 
excellent  review  are  available  (Hansen  and  Travis,  1974;  Shettle  and  Kuriyan, 
1975;  Lenoble,  1977).  This  section  will  focus  on  approximate  multiple 
scattering  methods  which  can  be  used  as  parameter izations  applicable  to  both 
LOWTRAN  and  FASCODE  models. 

A.l.l  Solar  Radiation 

For  practical  purposes  it  is  useful  to  categorize  the  extensive  hier¬ 
archy  of  available  computational  procedures  based  on  both  the  relative  accu¬ 
racy  of  the  results  required  and  the  expenditure  of  computational  resources 
incurred.  A  subjective  evaluation  according  to  these  criteria  results  in  tour 
classes  of  method  based  on  resultant  accuracy  characterized  as  either  exact  or 
approximate  and  coraputat"  al  methodologies  denoted  as  either  numerical  or 
analytical.  Exact  methous  are  defined  as  those  which  agree  among  themselves 
to  some  acceptable  degree  of  accuracy  when  applied  to  a  variety  of  diverse 
problems.  All  other  methods  are  approximate.  The  most  recent  intercomparison 
of  exact  methods  suggests  that  about  1%  accuracy  in  radiance  can  be  expected 
(Lenoble,  1977).  A  variety  of  exact,  numerical  methods  are  listed  in  Table 
A-l  along  with  representative  literature  citations.  The  essential  difference 
between  numerical  and  analytical  approaches  is  that  the  latter  generally 
require  less  computer  time.  Analytical  approaches  are  those  which  are 
generally  implementable  using  simple  algebraic  icldLions  supported  by  a  priori 
specifiable  parameters,  i.e.,  there  are  no  recurrent  or  iterative  calculations 
and  no  convergence  criteria.  Exact  analytical  methods  are  available  only  for 


a  few  rases  which  are,  unfortunately,  not  immediately  applicable  to  general 
terrestrial  radiance  simulation  situations.  These  cases  include,  for  the  most 
part,  approaches  based  on  Chandrasekhar's  H  functions  for  semi-infinite  atmo¬ 
spheres  (t  =  »)  and  X  and  Y  functions  for  finite  atmospheres  (Chandrasekhar, 
1960).  For  the  former  set  of  problems,  for  example,  solutions  are  available 
for  isotropic,  Rayleigh  and  various  anisotropic  phase  functions.  For  finite 
atmospheres,  solutions  are  available  only  for  isotropic  (Carlstedt  and  Mulli- 
kan,  1966)  and  Rayleigh  scattering  (Coulson  et  al.,  1960).  These  exact  treat¬ 
ments,  often  obtainable  in  the  form  of  lookup  tables,  are  useful  as  checks  for 
the  accuracy  of  both  exact  numerical  methods  and  proposed  approximate  methods 
when  applied  to  degenerate  cases  (such  as  isotropic  scattering). 

Approximate  methods  include  those  based  on  evaluating  the  first  few 
tractable  terms  or  orders  of  more  extensive  exact  numerical  treatments  and 
those  formulated  specifically  as  approximations.  In  the  former  class  of 
approximate  methods,  the  philosophy  is  to  accept  the  known  risk  of  increased 
error  for  a  specific  application  in  order  to  benefit  from  reduced  computer  run 
time.  Examples  include  terminating  successive  order  of  scattering  approaches 
after  a  few  scatterings  (see  Nagel  et  al.,  1978),  reducing  the  number  of 
Legendre  polynomials  carried  in  the  expansion  of  an  angular  scattering  func¬ 
tion,  and  limiting  the  number  of  quadrature  angles  in  the  discrete  ordinate 
method.  The  method  is  approximate  since  the  user  knows  a  priori  that  the  full 
available  power  of  the  exact  numerical  method  at  hand  is  not  being  exer¬ 
cised.  For  the  simplest  cases,  these  approaches  often  result  in  analytical 
solutions  which  may  themselves  be  useful  such  as  single  and  double  scattering 
solutions  (Deirraendj ian,  1969;  Hovenier,  1971),  the  two-  and  four-stream 
approximations  based  on  the  discrete  ordinate  method  (Liou,  1973,  1974),  and 
the  two-step  function  approach  based  on  general  variational  methods  (Sze, 

1976;  Burke  and  Sze,  1977). 

Multiple  scattering  methods  formulated  as  Intrinsically  approximate  have 
a  common  heritage  in  the  estimation  of  radiation  fluxes  within  larger 
numerical  models  (see  Stephens,  1984).  This  applies  both  to  the  use  of 
similarity  relations  to  reduce  anisotropic  problems  to  isotropic  equivalents 
(van  de  Hulst  and  Grossman,  1968)  and  various  so-called  two-stream 
approximations  (Meader  and  Weaver,  1980).  Notably,  it  is  only  recently  that 
various  investigators  have  sought  approaches  to  treat  the  approximation  of 
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multiply  scattered  radiances  by  augmenting  a  single  scattering  solution  with 
an  estimate  of  the  multiply  scattered  radiance  field  based  on  a  suitable  flux 
parameterization.  Included  among  these  are:  (.a)  the  "diffuse  field"  approx¬ 
imation  described  by  Bergstrom  et  al.  (1981)  for  use  in  visibility  impact  mod¬ 
els,  (b)  a  path  radiance  approximation  based  on  the  delta-Eddington  approach 
by  Hering  (1981)  which  was  implemented  within  an  operational  visible  contrast 
transmittance  model,  and  (c)  Kaufman's  (1979)  radiance  approximation  for  ap¬ 
plications  to  remote  sensing.  Kaufman  has  primarily  applied  his  approximate 
scheme  as  a  device  to  evaluate  multiply  scattered  contributions  to  upward 
radiance  within  the  context  of  multidimensional  remote  sensing  (cf.  Kaufman, 
1982).  Analogous  solutions  for  downward  radiance  were  derived  by  Isaacs 
(1981)  for  application  to  visibility  model  related  radiance  evaluation. 

This  overview  of  computational  methods  applicable  to  the  simulation  of 
radiances  due  to  the  multiple  scattering  of  solar  radiation  suggests  that  a 
suitable  parameterization  may  be  found  among  one  of  the  approximate  approaches 
described  above.  In  order  to  quantify  this  assessment,  radiance  comparisons 
were  made  between  representative  methods  from  each  of  the  categories  described 
above  for  a  set  of  model  atmospheres  with  aerosol  content  varying  from  ultra- 
clear  to  hazy  conditions.  Table  A-2  summarizes  the  scattering  treatments 
selected  in  this  comparison.  The  Gauss-Seidel  multiple  scattering  (MS)  ap¬ 
proach  was  selected  to  provide  exact  numerical  radiance  results  for  reference 
purposes  (Dave  and  Gazdag,  1970).  The  algorithm  used  (Dave,  1972)  evaluates 
fluxes  and  radiances  for  atmospheres  with  arbitrary  vertical  distribution  of 
anisotropic  aerosol  extinction  overlying  a  Lambert  reflecting  surface.  When 
applied  to  pure  molecular  scattering  cases,  radiance  results  from  the  code 
reproduce  the  exact  analytical  tabulations  of  Coulson  et  al.  (1960)  based  on 
Chandrasekhar's  X  and  Y  functions  for  Rayleigh  scattering  to  high  accuracy. 

Three  parameterization  approaches  were  investigated.  These  Include:  (a) 
single  scattering  (SS),  (b)  a  stream  approximation  (SA),  and  (c)  successive 
orders  of  scattering  (SOS).  The  SOS  method  (which  is  presently  solved  numer¬ 
ically)  was  formulated  to  treat  the  single  scattered  radiance  exactly  and  to 
approximate  higher  order  scattering  by  an  isotropic  angular  scattering  func¬ 
tion.  These  approaches  are  summarized  in  Table  A-3.  Results  of  the  test  case 
comparisons  are  discussed  in  Section  A. 2. 


Table  A-2 


Methods  Selected  for  Inter comparison 


Computational  Effort 


Accuracy  Numerical  Analytical 


Exact 

Gauss-Seidel  (MS] 

X,Y  functions* 

(Dave  and  Ga2dag,  1970; 

(Chandrasekhar,  1960; 

-  _  .. 

’Dave,  1972) 

Coulson  et  al.,  1960) 

Approximate 

Successive  Order 

Single  Scattering  [SS 

of  Scattering  [SOS] 

Stream  Approximation  fSA] 
(Kaufman,  1979) 


Source'  Function,  iJ,  anil  Radiance,  %  Evaluation  Approaches  for  Thr^'  Candidate  Approximation 
(8  ■  temperature,  T  -  transmittance,  t  =  optical  depth) 


A. 1.2  Thermal  Radiation 


The  detailed  numerical  schemes  described  in  the  previous  section  can  be 
used  to  perform  thermal  radiation  calculations  by  including  the  thermal 
emission  term.  For  example,  the  doubling  method  has  been  extended  by  Wiscombe 
to  include  an  inhomogeneous  internal  source  (see  Wiscombe,  1976  and  the  cited 
references  for  other  detailed  methods).  Although  these  methods  are  useful  for 
reference  calculations,  from  practical  considerations,  the  study  of  particle 
multiple  scattering  on  atmospheric  transmittance  has  to  rely  heavily  on 
■Oarameterlzatlons.  This  is  especially  true  for  the  LOWTRAN/FASCODE  models. 

Studies  of  the  effects  of  aerosol  multiple  scattering  on  atmospheric  in¬ 
frared  transmittance  have  emphasized  the  calculation  of  fluxes  rather  than  ra¬ 
diances.  The  two-stream  approach  has  been  the  commonlv-used  parameterization. 
In  general,  the  two-stream  approximation  is  a  method  of  representing  the  radi¬ 
ation  field  in  two  streams  in  one-dimensional  radiation  transfer  problems. 

The  approximation,  proposed  by  Schuster  (1905)  and  Schwarzschild  (1906),  is  to 
assume  that  the  intensity  in  the  positive  direction  is  isotropic  and  that  In 
-the  negative  direction  has  a  different  value  but  is  also  isotropic.  Chandra¬ 
sekhar  (1960)  has  used  the  Gaussian  quadrature  to  approximate  the  Integral 
over  angle,  terming  this  the  "discrete  ordinate"  method.  The  first  approxima¬ 
tion  of  the  discrete  ordinate  method  results  in  the  same  form  as  the  Rchuster- 
. Schwarzschild  approximation  except  for  different  constant  coefficients  due  to 
Gaussian  division. 

Sagan  and  Pollack  (1967)  and  Piotrowski  (1956)  compared  the  anproximat Ion 
with  the  exact  solutions  for  the  cases  of  conservative  and  non-conservative 
Isotropic  scattering.  The  results  show  5-10%  accuracy  in  absolute  values.  A 
survey  paper  by  Irvine  (1960)  has  investigated  several  approximation  methods, 
including  two-stream,  Eddington  and  modifed  two-stream  approximation.  Uses  of 
two-stream  approximation  for  treating  multiple  scattering  in  thermal  radiation 
are  numerous  (e.g.,  Rasool  and  Schneider,  1971;  Wang  and  Domoto,  1974;  Hansen 
et  al.,  1984),  This  parameterization  has  been  used  for  climate  application 
such  as  investigation  of  climatic  effect  due  to  aerosols  (Hansen  et  al.,  1980). 


A . 2  Trade-off  Analysis  of  Exact  an d  Approximate  Treatments 


A  series  of  computer  codes  was  written  to  support  an  intercomparison  of 
angle  dependent  radiances  from  single  scattering  (SS),  stream  approximation 


(SA),  and  successive  orders  of  scattering  (SOS)  parameterizatl ons  with  those 
evaluated  using  accurate  numerical  multiple  scattering  treatments.  For  the 
stream  approximation,  two  independent  flux  treatments  were  employed,  one  for 
solar  radiation  and  one  for  thermal  radiation.  The  intercomparison 
calculations  were  performed  by  parametrically  varying  the  relevant  optical 
properties  of  optical  depth,  single  scattering  albedo,  and  asymmetry  factor. 
Thus,  they  were  offline  from  either  the  LOWTRAN  or  FASCODE  models. 
Furthermore,  these  calculations  assumed  single  homogeneous  scattering  layers. 

Assessment  criteria  include  accuracy  defined  in  terms  of  computer 
execution  times  for  a  set  of  test  cases  and  accuracy  measured  by 
intercomparison  with  numerical  results. 


A , 2 . 1  Solar  Scattering 

For  solar  scattering,  six  model  atmospheres  were  used  in  a  comparison 
test  consisting  of  one  Rayleigh  scattering  case  and  five  atmospheres  including 
increasing  amounts  of  aerosol.  The  vertical  distribution  of  molecular  density 
was  that  of  the  AFGL  midlatitude  summer  model  while  that  for  the  aerosol  was 
Elterman's  (1970)  distribution.  The  turbidity  weighting  factor  (defined  as 
the  ratio  of  total  aerosol  scattering  optical  depth  to  total  scattering 
optical  depth)  varied  from  zero  to  near  unitv  with  corresponding  surface 
visual  ranges  decreasing  to  about  1.2  km  in  haze.  The  optical  properties  of 
these  model  atmospheres  are  given  in  Table  A-4.  For  each  of  the  six  test 
cases  radiances  were  calculated  using  the  SS,  SOS,  and  SA  methods  at  nine 
values  of  zenith  angle,  three  values  of  azimuth  angle,  and  two  surface  reflec¬ 
tances  for  both  upward  viewing  from  the  surface  and  downward  viewing  from 
space  (i.e.,  a  total  of  108  values).  The  solar  zenith  angle  for  all  cases  was 
60®.  Corresponding  exact  multiple  scattering  (MS)  cases  were  evaluated  using 
the  Gauss-Seidel  iterative  method  (Dave,  1972).  The  phase  function  used  was 
based  on  the  radiance  values  from  the  exact  method,  R(MS)  and  each  approximate 
method  [R(A),  A-SS,  SOS,  SA],  percent  errors  were  calculated  for  the 
corresponding  ith  viewing  angle,  surface  reflectance  pair  according  to: 


P  F.  1  ( A ) 


r R(MS)-R( A) 
1  K  ( M  S ) 


X  100"/ 


(A.  I  ) 


Additionally  in  order  to  provide  a  measure  of  overall  accuracy  for  each  case, 
root  mean  square  errors,  RMS,  were  evaluated  for  each  approximate  method.  A, 
and  surface  reflectance,  r,  according  to: 

54  1/2 

RMS(A,r )  -  (  I  (PE  (A) I  )  (A. 2) 

i-1 

A  comparison  of  RMS  error  according  to  Equ.  (A. 2)  is  summarized  in  Table 
A-5,  For  each  of  the  approximations  (single  scattering,  stream  approximation, 
successive  orders  of  scattering),  RMS  errors  are  provided  for  the  six  model 
atmospheres  and  two  selected  values  of  surface  reflectance  (r  *  0,0,  0.4), 

The  magnitude  of  RMS  errors  reported  in  Table  A-5  for  single  scattering 
approximate  those  which  would  be  obtained  using  the  LOWTRAN6  single  scattering 
option  (IEMSCT  *  2)  for  a  like  sun/observer  geometry.  Examining  Table  A-5,  it 
is  notable  that  the  increase  of  RMS  error  with  optical  thickness  is 

* 

considerably  less  for  the  stream  approximation  and  successive  scattering 
method  than  for  the  single  scatter  parameterization.  This  is  not  surprising 
since  the  latter  does  not  treat  the  multiply  scattered  radiance  contribution 
which  becomes  increasingly  Important  at  larger  optical  thicknesses.  At  a 
moderate  optical  thickness  (say  1.0,  corresponding  to  an  approximate  visual 
range  of  6.4  km)  and  zero  surface  reflectance  the  SOS  and  SA  approaches  are 
better  by  factors  of  about  1.5  and  2.0,  respectively. 

In  order  to  provide  a  measure  of  efficiency  for  each  approximation,  exe¬ 
cution  times  were  tabulated  for  calculation  of  scattered  radiances  for  the 
ensemble  of  eighteen  upward  and  downward  zenith  angles,  three  azimuth  angles, 
and  two  values  of  surface  reflectance.  These  results  are  presented  In  Table 
A-6.  The  codes  to  evaluate  each  approximation  were  run  on  AER's  Harris  H800 
and  while  not  optimized  for  speed  were  configured  for  consistency  with  the 
layered  atmosphere  code  structure  of  the  L0WTRAN/FASC0DE  models.  The  exact 
multiple  scattering  (MS)  execution  times  were  obtained  from  comparison  runs 
mad  •  on  the  APGL  CYBEK.  To  provide  a  basis  for  comparison,  Harris  equivalent 
execution  times  (the  values  in  parentheses)  were  estimated  for  each  case  using 
a  benchmark  execution  of  the  Dave  code  from  each  machine  to  obtain  a  conver¬ 
sion  factor  of  about  2.1.  It  is  Important  to  note  that  these  times  apply  to 
the  evaluation  of  the  radiances  from  prespecified  profiles  of  optical  param¬ 
eters  and  do  not  include  time  spent  calculating  atmospheric  transmission. 


Table  A-5 


Comparison  RMS  Errors  for  Single  Scattering  (SS), 
Stream  Approximation  (SA) ,  and  Successive 
Orders  of  Scattering  (SOS) 


Optical  Thickness/ 
Ref lectivlty 


0.10 

(Rayleigh) 

/0.0 

/0.4 

0.25 

(Case  l) 

/0.0 
/0 .4 

0.50 

(Case  2) 

/0.0 
/0 .4 

1.00 

(Case  3) 

/0.0 
/0 .4 

2.00 

(Case  A) 

/0.0 
/0 .4 

5.00 

(Case  5) 

/0.0 
/0 .4 

SS 

SA 

SOS 

13.5 

4.3 

4.5 

29.2 

10.4 

3.6 

33.4 

12.2 

17.9 

41.8 

4.9 

18.8 

50.4 

17.2 

29.5 

57.6 

9.4 

31.9 

67.0 

29.1 

43.3 

72.5 

23.2 

48.6 

80.8 

34.0 

58.1 

83.3 

30.3 

63.5 

90*0 

39.1 

73.1 

90.5 

40.1 

73.7 

It  is  interesting  to  compare  execution  times  relative  to  the  correspond¬ 
ing  single  scattering  values.  The  stream  approximation  takes  about  two  times 
longer  regardless  of  optical  thickness.  This  is  of  interest  since  the  succes¬ 
sive  scattering  execution  times  increase  with  increasing  optical  thickness  and 
are  at  best  eight  times  larger,  although  this  could  probably  be  shortened  by 
careful  optimisation*  Except  for  the  Rayleigh  case  (which  could  be  solved 
exactly  anyway),  multiple  scattering  execution  times  are  at  least  an  order  of 
magnitude  larger  and  Increase  slightly  with  optical  thickness. 

A . 2 , 2  Thermal  Scattering 

A  set  of  model  atmosphere  test  cases  was  employed  to  compare  each  mul¬ 
tiple  scattering  parameterization  for  thermal  scattering  with  an  exact  numer¬ 
ical  treatment.  Spectral  regions  of  interest  included  the  infrared  where 
scattering  by  aerosol  and  cloud  can  be  a  factor  in  determining  the  radiance 
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Table  A-6 


Comparison  Execution  Times  for  Single  Scattering  (SS), 
Stream  Approximation  of  (SA),  and  Successive 
Order  of  Scattering  (SOS)  Methods 


Execution  Time  (sec) 


Optical  Thickness 

SS1 

SA1 

SOS1 

MS2'3 

0.10  (Rayleigh) 

3.4 

7.8 

27.4 

2.45 

(5.1) 

0.25  (Case  1) 

3.4 

7.8 

37.6 

47.64 

(100.1) 

0.50  (Case  2) 

3.4 

7.8 

47.6 

48.79 

(102.5) 

1.00  (Case  3) 

3.4 

7.8 

63.4 

50.14 

(105.3) 

2.00  (Case  4) 

3.4 

7.8 

90.6 

51.98 

(109.2) 

5.00  (Case  5) 

3.4 

7.8 

107.8 

52.92 

(111.13) 

Harris  H800  CPU  time  to  evaluate  radiances  for  eighteen  upward  and  downward 
zenith  angles,  three  azimuth  angles  and  two  values  of  surface  reflectance, 
i.e.,  a  total  of  108  atmospheric  paths. 

2AFGL  CYBER  CPU  execution  time  in  seconds  for  multiple  scattering  (Dave,  1976) 
solution  of  same  108  atmospheric  paths. 

^Parentheses  contain  estimated  Harris  H80Q  execution  times  for  comparison  (the 
Rayleigh  scattering  case  is  actual)  based  on  «n  H800/CYBER  time  factor  cf  2.1 
obtained  from  a  common  benchmark  run  of  the  Dave  code. 


field  and  the  microwave/mi 11 (me  ter  wave  where  scattering  is  due  largelv  to 
precipitation.  Calculations  were  performed  for  similar  sets  of  optical 
properties  for  the  stream  approximation  (SA),  successive  orders  of  scattering 
(SOS),  and  an  accurate  numerical  method  for  multiple  scattering,  the  discrete 
ordinate  method  (DOM). 

In  order  to  maximize  the  effects  of  aerosol  scattering  in  the  infrared,  a 
spectral  region  was  selected  (900  cm-1)  where  gaseous  absorption  is  relatively 
weak.  This  window  region,  influenced  primarily  by  water  vapor  continuum 
absorption,  is  used  for  remote  sensing  purposes  such  as  surface  temperature 
retrieval.  The  test  comparisons  were  done  parametrically  with  optical  thick¬ 
nesses  in  the  range  from  0.05  to  1.00,  single  scattering  albedos  between  0.5 
and  o,9,  and  asymmetry  factors  between  0.6  and  0.8.  These  ranges  were  chosen 
both  to  include  values  representative  of  the  10-12  urn  window  for  typical  haze 
conditions  and  to  push  the  effect  of  scattering  to  unrealistically  high  values 
for  testing  purposes.  The  optical  properties  of  the  model  atmospheres  used  in 
the  thermal  scattering  comparison  are  summarized  in  Table  A-7. 


Table  A-7 

Optical  Properties  of  Model  Atmospheres 
Used  in  Thermal  Scattering  Comparison 

Optical  Single  Scattering  Asvmmetrv 

Thickness  Albedo  Factor 


Infrared  Region  (v  *  900 

cm”1  ) 

Case  1  a,b 

0.05 

0.50, 

0.75, 

C.90 

0.6,  0.7,  0.8 

Case  2  a,b 

0.25 

0.50, 

0.75, 

0.90 

0.6,  0.7,  0.8 

Case  3  a,b 

0.50 

0.50, 

0.75, 

0.90 

0.6,  0.7,  0.8 

Case  4  a,b 

1.00 

0.50, 

0.75, 

0.90 

0.6,  0.7,  0.8 

Microwave/millimeter  wave  (f  = 

19.35,  37.00  GHz) 

Case  5  a,b  (19.35  GHz) 

0.906 

(M-P,  Rain 

0.21 

rate 

15  mmhr"1) 

Phase  Function 

Case  6  a,b  (37.00  GHz) 

3.20 

(M-P,  Rain 

0.39 

rate 

1 5  mmh  r - * ) 

Phase  Function 

The  surface  temperature  for  all  calculations  was  assumed  to  be  288  K„ 

Two  types  of  atmospheric  temperature  profile  were  Investigated:  an  Isothermal 
atmosphere  with  a  temperature  of  283  K  and  an  atmosphere  with  a  10  K  tempera¬ 
ture  gradient  from  top  to  bottom  and  a  mean  temperature  of  283  K.  These  are 
called  the  isothermal  and  gradient  atmospheres,  respectively. 

Microwave/millimeter  wave  calculations  were  done  for  cases  with  similar 
surface  temperatures  and  temperature  profiles.  Two  frequencies,  19.35  and 
37.0  GHz  (1.55  and  0.81  cm,  respectively)  corresponding  to  those  used  for 
satellite  remote  sensing  were  investigated.  For  these  cases  (see  Table  A-7), 
scattering  was  assumed  due  to  a  Marshall-Palmer  distribution  of  rain  with  a 
rain  rate  of  15  mm/h.  Precipitation  scattering  properties  such  as  the  extinc¬ 
tion  coefficient,  single  scattering  albedo,  and  angular  scattering  fuction  were 
evaluated  from  a  parameterization  of  the  Mie  theory  calculations  of  Savage 
(1978)  developed  for  a  millimeter  wave  simulation  model  (Isaacs  et  al.,  1985). 

The  discrete  ordinate  method  (DOM)  (Liou,  1973)  provided  numerical 
solutions  as  a  standard  for  comparison,  A  sixteen  stream  version  of  the  DOM 
code  was  used.  Radiance  results  were  thus  given  at  eight  upward  and  eight 
downward  zenith  angles  within  the  domains  from  0  to  90  degrees  (upward 
radiance)  and  90  to  180  degrees  (downward  radiance).  The  endpoints  of  these 
ranges  (i.e.,  0,  90,  and  180  degrees)  are  extrapolated  from  internal  points 
and  are  thus  not  accurate  values.  To  facilitate  calculation  of  the  required 
Legendre  polynomial  expansion  of  the  angular  scattering  function,  Henyey- 
Greenstein  phase  functions  with  the  desired  asymmetry  factors  were  employed. 

Results  of  the  multiple  scattering  parameterization  intercomparison  were 
tabulated  as  radiances,  percent  errors  as  a  function  of  zenith  angle,  and  per¬ 
cent  root  mean  squared  (RMS)  errors.  The  latter  quantity  was  evaluated  over 
all  zenith  angles  used  in  the  intercomparison.  Sample  results  are  shown  in 
Figures  A-l  to  A-4.  Figure  A-l(a)  and  (b)  illustrate  the  comparison  between 
upward  and  downward  emergent  radiances,  respectively,  for  the  SA  and  SOS  par- 
ameterization8  and  the  exact  DOM  solution.  This  particular  calculation  is  for 
an  optical  depth  of  0.25,  single  scattering  albedo  of  0.9,  and  an  asymmetry 
factor  of  0.8.  The  SOS  approach  underestimates  upward  radiance  and  overesti¬ 
mates  downward  radiance  compared  to  the  DOM  due  to  the  assumption  of  isotropic 
scattering.  The  SA,  however,  reasonably  reproduces  the  multiply  scattered 
radiance  field,  except  near  the  horizon  (zenith  angles  approaching  90°). 
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Percent  RMS  erros  for  SA  and  SOS  compared  to  DOM  as  a  function 
of  optical  depth  and  single  scattering  albedo  for  g  *  0,6: 

(a)  SA  upward  radiance,  (b)  SOS  upward  radiance,  (c)  SA  downward 
radiance,  and  (d)  SOS  downward  radiance. 
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Figure  A-4,  Percent  RMS  erros  for  SA  and  SOS  compared  to  DOM  as  a  function 
of  optical  depth  and  single  scattering  albedo  for  g  -  0.8: 

(a)  SA  upward  radiance,  (b)  SOS  upward  radiance,  (c)  SA  downward 
radiance,  and  (d)  SOS  downward  radiance. 


•—Similar  conclusions  can  be  drawn  by  looking  at  the  internal  upward  and 

-V  ,¥  ■ 

--downward  radiances  halfway  through  the  layer  at  an  optical  depth  of  0.125 
(Figures  A-2(a)  and  (b),  respectively). 

-  'ffr  Figures  A-3  and  A-4  illustrate  the  percent  RMS  error  for  the  stream  ap¬ 
proximation  and  successive  orders  of  scattering  compared  to  the  discrete  ordi- 
••  .nate  solutions  for  all  optical  depths,  and  single  scattering  albedos,  and 
v*5»scattering  asymmetry  factors  of  0.6  and  0.8,  respectively.  While  the  SOS 
ijagfiapproach  provides  reasonable  upward  radiances  with  errors  only  slightly  great- 
than  the  SA,  SOS  derived  downward  radiances  are  quite  inaccurate.  Downward 
_  :’radiance8  are  a  real  test  of  the  scattering  parameterization  since  upward  ra- 
'  diances  are  largely  determined  by  the  uniform  surface  emission  at  these  optic¬ 
al  depths.  Notably,  RMS  errors  for  the  stream  approximation  are  weakly  depen¬ 
dent  on  optical  depth  and  are  at  most  about  25  percent.  Errors  are  largest 
?3for  high  single  scattering  albedo  and  asymmetry  factor,  i.e.,  when  the  scat- 
altering  effect  is  maximum  and  when  the  angular  scattering  function  is  strongly 
1  forward  scattering.  This  is  expected  due  to  the  nature  of  the  approximation 
raised  • 

•  •ip  "•**■  - 

Computer  execution  times  were  also  tabulated  for  each  thermal  scattering 
test  case  ensemble  to  ascertain  their  relative  efficiency.  Execution  times 
iare  given  in  Table  A-8.  Times  for  the  SOS  and  DOM  are  comparable  while  the  SA 
^iwas  about  40  times  faster  than  the  DOM.  These  results  suggest  that  the  SA  is 
•-.t.-.a^he  candidate  parameterization  of  choice.  . r  ' 


■  :  ?  A. 3  Recommendations  for  LOWTRAN  and  FASCODE 

'g  m  i  _  tu  -  i  ■_  ■  m 

;  ■  -  Based  on  these  results,  the  stream  approximation  was  selected  for 
■' implementation  for  thermal  scattering  as  well  as  solar  scattering  in  the 
; 0  WT RAN  and  FASCODE  models.  Based  on  the  evaluation  and  consideration  of  the 
requirements  of  efficiency,  accuracy  and  code  structure  of  LOWTRAN  and  FASCODE 
models,  the  parameterization  of  the  stream  approximation  together  with  the 
■  -  correlated-k  approximation  is  best  suited  for  implementation  into  these 
'"models.  The  parameterization  is  unique  in  the  sense  that  the  correlated-k 
approximation  decouples  the  multiple  scattering  from  the  treatment  of  nongray 


Table  A-8 

Therm  1  Results 


Execution  Time  (sec.) 
Optical  Thickness 


SA 


SOS 


DOM  (Liou)* 


_  ,/?Tau  *  0.05 


<g  “  0.6,  o)q  ■  0.50) 

|  (12.10) 

-  0.6,  a>0  -  0.90) 


.-.r,T.  --.Tau  -  0.25 

-  0.6,  u0  -  0.50) 

(g  *  0.6,  (o q  *  0.90) 

— r^^ESTau  -0.50 

•  : “  0.6,  (i)q  *  0.50) 

~'~~?rr  *'-~(g  *  0.6,  uq  ”  0.90) 


'  5^  ■* 1  ■  ” 


•"V  -  • 


Tau  -  1.00 

‘  (g  *  0.6,  wq  «  0.50) 
(g  ■  0.6,  Wf.  *  0.90) 


f;fisr  -V. 


i-  .35 

6.28 

5.76 

f  .35 
(12.10) 

^6.28 

5.79 

1  .34 

:<  1 2 . 1 6 ) 

7.55 

5.79 

.34 

.(12.19) 

T 

10.17 

5,81 

~f-r 

.34 

(12.02) 

8.87 

5.72 

.34 

(12.15) 

12.79 

5.79 

.34 

(12.06) 

10.14 

5.74 

.36 

(12.22) 

16.69 

5.82 

.V' 


\%Runs  were  made  on  AFGL’s  Cyber  750  system.  All  other  runs  on  AER's  Harris 
5^1800  system.  Times  In  parenthesis  are  Harris  equivalents  for  runs  made  on 
t-The  Cyber,  using  a  benchmark  factor  of  2.1. 


gaseous  absorption  so  that  the  stream  approximation  can  be  used  to  treat 
multiple  scattering  in  both  the  LOWTRAN  and  FASC0DF.  models.  The  stream 
approximation  is  extremely  attractive  since  it  retains  the  analytical 
simplicity  (and  hence,  efficiency)  of  single  scattering  while  improving  the 
treatment  of  multiple-scattering  (and  hence,  accuracy). 
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Appendix  B 

EXAMPLE  OF  FASCODE  OUTPUT 

A  case  is  Included  as  a  test  for  FASCODE  with  multiple  scattering.  The 
input  cards  for  the  case  are  listed  in  Table  B-l  and  the  listing  of  the  file 
TAPE  6  is  in  Table  B-2.  The  next  section  will  discuss  the  output  for  the  test 
case.  Only  those  details  which  differ  from  the  standard  FASCODE  output  will 
be  discussed. 

This  case  demonstrates  the  new  multiple  scattering  calculation  for 
thermal  scattering  of  downward  radiance.  The  parameters  selected  for  this 
ji ase  are  as  follows:  the  atmospheric  profile  is  the  1962  U.S.  Standard  and 
the  boundary  layer  aerosol  model  is  RURAL  with  5  km  meteorological  range.  The 
path  is  a  slant  path  from  the  ground  to  20  km  with  a  zenith  angle  at  the 
ground  of  0.0  degrees. 

The  output  for  this  case  shown  in  Table  B-2  will  now  be  described.  The 
first  difference  is  in  CARD  2,  which  now  contains  another  quantity,  IMS.  IMS 
Sis  the  eleventh  variable  printed  on  CARD  2,  and  is  the  flag  for  multiple 
scattering.  When  IMS  is  1,  the  multiple  scattering  features  of  the  program 
will  be  used.  The  next  difference  which  occurs  in  the  output,  is  a  print  of 
the  path  parameters  which  occurs  on  pages  2  and  3.  This  occurs  due  to  a  call 
to  FSCGEO  which  is  needed  to  determine  whether  or  not  the  user  defined  path 
starts  at  or  intersects  the  ground,  and  whether  or  not  HMINMS  and  HMAXMS  are 
contained  within  the  user  defined  path.  (This  is  explained  further  in  the 
user  instructions  in  Section  3,1.4.  If  additional  paths  were  determined  to  be 
necessary  in  order  to  accommodate  the  multiple  scattering  calculation,  these 
paths  would  be  run  and  the  corresponding  output  would  appear  in  the 
calculations  at  this  point.  The  user  defined  path  would  be  the  last  path  that 
would  appear  in  the  output.) 

The  next  25  pages  in  the  output  shown  in  Table  B-2  are  identical  to  those 
in  the  output  for  a  run  without  multiple  scattering.  The  following  21  pages 
which  output  the  optical  depths,  radiance,  and  transmittance  for  each  FASCODE 
layer  are  different  for  the  multiple  scattering  runs.  The  differences  will  be 
discussed  for  one  layer,  since  each  layer  has  the  same  output.  The  first 
difference  is  that  the  radiance  that  is  printed  is  the  thermal  emission  and 
not  the  conservative  scattering  radiance  supplied  with  the  standard  FASCODE 


-run.  The  second  difference  is  that  after  the  radiance  and  transmittance  is 
printed,  the  composite  upward  flux  and  reflectance  for  that  layer  is  also 
printed.  These  quantities  are  intermediate  results  used  by  the  multiple 
scattering  sections  of  the  code,  and  are  printed  for  the  users  benefit  to  see 
what  the  multiple  scattering  is  doing. 

That  would  conclude  the  output  for  the  standard  FA8C0DE  run,  but  for  the 
multiple  scattering  run,  we  still  need  to  add  the  fluxes  from  the  top  down  in 
order  to  determine  the  multiple  scattered  contribution.  The  next  21  pages  are 
printout  from  the  downward  adding  of  the  fluxes.  The  first  page  starts  off 
—with  the  composite,  downward  flux  and  reflectance  for  layer  10.  Then  the  total 
upward  and  downward  flux  for  the  top  of  layer  10  is  printed.  The  next  page 
prints  out  the  composite  downward  flux  and  reflectance  for  layer  9.  Then  the 
total  upward  and  downward  flux  for  the  top  of  layer  9  is  printed.  The  next 
page  starts  with  the  output  of  the  source  function  for  the  multiple  scattering 
contribution  to  the  radiance  from  Layer  10  (this  quantity  is  actually  the  MS 
=semi8Bion  from  layer  10),  This  is  followed  by  the  composite  downward  flux  and 
reflectance  for  layer  8,  followed  by  the  total  upward  and  downward  flux  for 
the  top  of  layer  8.  The  source  function  for  layer  9  appears  at  the  top  of  the 
next  page,  and  is  followed  by  the  total  multiple  scattered  emission  from  both 
layer  9  and  10.  This  output  scheme  Is  continued  such  that  on  the  last  of  the 
21  pages  the  final  MS  emission  that  Is  printed  Is  the  total  multiple 
scattering  emission  for  the  path. 

The  next  and  final  page  In  the  output  is  the  total  radiance  for  the 
FASCODE  run  and  consists  of  the  multiple  scattered  emission  added  to  the 
thermal  emission. 
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Appendix  C 

EXAMPLES  OF  LOWTRAN  OUTPUT 


Four  cases  are  included  as  tests  for  LOWTRAN  with  multiple  scattering. 
The  input  cards  for  the  four  cases  are  listed  in  Table  C-l  and  the  listing  of 
the  file  OUTPUT  (TAPE  6)  is  in  Tables  C-2  -  C-5,  The  next  section  will 
discuss  the  output  for  case  1.  Only  those  details  which  differ  from  the 
standard  LOWTRAN  output  will  be  discussed. 

C. 1  Case  1:  Thermal  Scattering  of  Downward  Radiance 

This  case  demonstrates  the  new  multiple  scattering  calculation  for 
thermal  radiance.  The  parameters  selected  for  this  case  are  as  follows:  the 
atmospheric  profile  is  the  1962  U.S.  Standard  and  the  boundary  layer  aerosol 
model  is  RURAL  with  5  km  meteorological  range.  The  path  is  a  slant  path  from 
the  ground  to  20  km  with  a  zenith  angle  at  the  ground  of  0.0  degrees. 

The  output  for  this  case  shown  in  Table  C-2  will  now  be  described.  The 
first  difference  is  in  CARD  1,  which  now  contains  another  quantity,  1MULT. 
1MULT  is  the  last  variable  printed  on  CARD  l,  and  is  the  flag  for  multiple 
scattering.  When  I  MULT  is  1,  the  multiple  scattering  features  of  the  program 
will  be  used.  A  message  to  this  effect  is  printed  in  the  output  directly 
after  the  message  which  tells  the  user  that  'PROGRAM  WILL  COMPUTE  RADIANCE'. 

The  next  four  pages  which  print  out  the  'ATMOSPHERIC  PROFILES'  and  the 
'REFRACTED  PATH  THROUGH  THE  ATMOSPHERE'  are  Identical  to  those  printed  out  by 
LOWTRANo.  The  following  three  pages  are  not  included  in  L0WTRAN6 .  These 
three  pages  contain  the  new  slant  path  parameters  for  using  multiple 
scattering  followed  by  the  refracted  path  used  by  the  multiple  scattering  run 
and  the  nummary  of  the  multiple  scattering  geometry  calculation.  The  final 
page  is  exactly  Identical  except  that  the  radiance  that  is  printed  out  is  the 
multiple  scattered  radiance. 

C.2  C  a  u c  2:^  Thermal  Scattering  ol  Upward  Had  lance 

This  cnse  Is  directly  analogous  to  Case  1  except  that  the  observer  Is 
looking  straight  down  instead  of  straight  up.  The  output  differences  are  tie 
same  -is  those  <_xp  1  ••  I  ried  lot  Cose  1. 


C . 3  Case  3:  Solar  Scattering  of  Downward  Radiance 


Case  3  illustrates  the  multiple  scattered  solar  radiance  calculation. 

The  parameters  for  this  case  are  the  same  as  those  for  Case  1,  with  the 
addition  of  the  sun  at  60  degrees.  The  output  differences  are  the  same  as 
those  listed  for  Case  1,  excluding  the  final  radiance  output. 

The  final  radiance  output  has  been  expanded  for  the  solar  case.  Two  new 
columns  have  been  added.  Under  the  column  heading  'PATH  SCATTERED  RADIANCE', 
there  are  now  three  columns.  These  are  'TOTAL'  (in  both  cm"1  and  microns)  and 
'S  SCAT'  (in  cm-1).  TOTAL  is  the  total  path  scattered  radiance,  containing 
both  single  and  multiple  scattered  contributions.  S  SCAT  is  the  single 
scattered  part  of  the  total  path  scattered  radiance.  This  S  SCAT  is  the 
quantity  that  LOWTRAN  would  provide  under  the  'PATH  SCATTERED'  heading  if  the 
user  were  to  run  LOWTRAN  without  multiple  scattering.  Likewise  under  the 
header  'GROUND  REFLECTED  RADIANCE'  there  are  now  three  columns.  These  are 
TOTAL  (in  both  cm"1  and  microns)  and  'DIRECT'  (in  cm"1).  TOTAL  is  the  total 
ground  reflected  contribution  to  the  radiance,  containing  both  scattered  and 
direct  contributions,  DIRECT  is  the  direct  ground  reflected  part  of  the  total 
reflected  radiance.  The  DIRECT  term  is  the  quantity  that  LOWTRAN  would 
provide  under  the  GROUND  REFLECTED  heading  if  the  user  were  to  run  LOWTRAN 
without  multiple  scattering.  The  significance  of  the  S  SCAT  and  the  DIRECT 
terms  is  that  the  user  can  run  one  case  of  multiple  scattering  and  by 
comparing  these  quantities  to  the  total  quantities  can  determine  the  relative 
importance  of  multiple  scattering  for  that  case.  This  will  allow  the  user  to 
decide  more  easily  whether  or  not  multiple  scattering  has  a  significant 
contribution  for  that  particular  case. 

C.4  Case  4:  Solar  Scattering  of  Upward  Radiance 

This  case  Is  directly  analogous  to  Case  3  except  that  the  observer  Is 
looking  straight  down  instead  of  straight  up.  The  output  differences  are  the 
those  expl  nl  lied  for  Case  3. 
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Test  Case  1 


INPUT  CARDS  FOR  THE  FOUR  LOWTRAN  TEST  CASES 
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